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A model of an elastic manifold driven through a random 
medium by an applied force F is introduced and studied. The 
focus is on the effects of inertia and elastic waves, in partic- 
ular stress overshoots in which motion of one segment of the 
manifold causes a temporary stress on its neighboring seg- 
ments in addition to the static stress. Such stress overshoots 
decrease the critical force for depinning and make the depin- 
ning transition hysteretic with static and pinned configura- 
tions coexisting with the steadily moving phase for a range 
of F. We find that the steady state velocity of the moving 
phase is, nevertheless, history independent and the critical 
behavior as the force is decreased is in the same universality 
class as in the absence of stress overshoots — the dissipative 
limit in which hysteresis cannot occur and theoretical anal- 
ysis has been possible. To reach this conclusion, finite-size 
scaling analyses have been performed and a variety of quan- 
tities studied, including velocities, roughnesses, distributions 
of critical forces, and universal amplitude ratios. 

If the force is increased slowly from zero, the behavior is 
complicated with a spectrum of avalanche sizes occurring that 
seems to be quite different from the dissipative limit. Related 
behavior is seen as the force is increased back up again to 
restart the motion of samples that have been stopped from 
the moving phase. The restarting process itself involves both 
fractal and bubble-like nucleation. Hysteresis loops in small 
and intermediate size samples can be understood in terms of 
a depletion layer caused by the stress overshoots. Surpris- 
ingly, in the limit of very large samples the hysteresis loops 
vanish. Although complicated crossovers complicate the anal- 
ysis, we argue that the underlying universality class govern- 
ing this pseudohysteresis and avalanches is again that of the 
apparently-very-different dissipative limit. But there are his- 
tory dependent amplitudes — associated with the depletion 
layer — that cause striking differences over wide ranges of 
length scales. Consequences of this picture for the statistics 
and dynamics of earthquakes on geological faults are briefly 
discussed. 



I. INTRODUCTION 

Extended elastic manifolds pulled through a quenched 
random medium by an applied force F exhibit, in the ab- 
sence of thermal fluctuations, a sharp transition from a 
pinned phase to a moving phase as F is increased through 
a critical value F c 0]. Examples include interfaces be- 
tween two fluids in porous media || or between oppo- 
sitely magnetized ferromagnetic domains, vortex lines 
and lattices in type II superconductors j|, [Q, charge 



density waves ^ , and planar crack fronts in solids || . 

Although the depinning transitions of interest are 
driven non-equilibrium transitions, it is instructive to 
draw an analogy with equilibrium phase transitions with 
the average velocity v playing the role of an order param- 
eter, F a tuning parameter, and the quenched variations 
of the random potential loosely analogous — at least as 
giving rise to an ensemble — to thermal fluctuations. 
The character of depinning transitions can, one might 
expect, be either discontinuous transitions with hystere- 
sis — loosely like first-order transitions — , or critical 
- analogous to second-order — transitions depending 
on the system and, perhaps, on its history; such history 
dependence is an effect that cannot occur in true equilib- 
rium. Theoretical analysis has shown that a broad class 
of realistic models undergo a critical depinning transition 
with a unique, history independent critical force in the 
limit of a large system and 

v ~ (F - F c f (1) 

for F just above F c . Several different universality classes 
have been studied, including both short and long-range 
interactions, and random forces with or without period- 
icity — the former arising for manifolds with a periodic 
structure in the direction in which they move || , |7j , || , 

But most of the theoretical analysis has focused on dis- 
sipative dynamics for which both inertia and any wave 
or other non-local stress propagation effects are ignored. 
The purpose of this paper is to study some of the conse- 
quences of these and other effects which we shall generi- 
cally refer to, for reasons to be explained shortly, as stress 
overshoot effects. 

In order to understand the potential of these effects to 
significantly change the nature of depinning, it is neces- 
sary to consider the nature of the irregular local motion 
that underlies the critical depinning phenomena. The 
elasticity of the manifold mediates between two compet- 
ing types of forces: the applied driving force and the 
random local pinning forces. For small F, the pinning 
dominates and the system relaxes to one of many static 
locally-stable configurations. But as the force is slowly 
increased, there will be local instabilities when the driv- 
ing force exceeds the random pinning in some small re- 
gion. A segment of the manifold will then move forward 
rapidly and there will be some transient motion, lim- 
ited in spatial extent, until a new static configuration is 
reached. For large F, in contrast, the applied force will 
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dominate and the system will approach a nonequilibrium 
statistically steady state with a non-zero mean velocity. 
Nevertheless, especially if the force is not too large, the 
motion on short length and time scales will be very ir- 
regular with instantaneous local velocities that far exceed 
the macroscopic average velocity v. 

In the absence of inertia or wave propagation, each seg- 
ment of the manifold will move in response to the total 
force applied to it: from the applied drive, from the ran- 
dom pinning, and from the other segments via elasticity. 
As long as the applied force is non-decreasing in time, 
the motion (at least after initial transients have decayed 
away) will be only in the "forward" direction in which 
the system is driven. This, combined with the convexity 
properties of elasticity, means that the static configura- 
tion the system will settle into after it is disturbed by an 
increase — global or local — in the applied force does not 
depend on the details of its dynamics. The motion can 
thus be considered as quasistatic, in spite of the rapid 
local motion that occurs. 

Now consider what can happen in the presence of either 
inertia or elastic waves that carry stress from one region 
to another in response to motion of one segment. The 
local dynamics would then appear to be crucial: If the 
local motion is rapid enough that the relaxation to a new 
static configuration is underdamped, a moving segment 
can overshoot one or more potential static configurations 
before settling, if at all, in another. Even if the inertia 
is small enough that such local prolonged jumps do not 
occur, as long as the motion is sometimes underdamped, 
a segment can temporarily overshoot a static configu- 
ration before relaxing back into it; this will produce a 
temporary overshoot of the stress — above its eventual 
static value — that this motion induces on neighboring 
segments. Any arbitrary small overshoot in the stress 
has the potential to dislodge another segment if there is 
one nearby that was sufficiently close to being destabi- 
lized in the absence of the overshoot; again, the effects 
of this will be to cause the system to skip through a po- 
tentially static configuration without stopping. Elastic 
waves, just like their electrodynamic cousins, carry with 
them pulses in stress that are larger than the eventual 
static stress that will obtain long after the waves have 
passed by. These stress overshoots, like those from the 
inertia of local motion, have the potential to cause over- 
jumping. 

Very generally, overjumping of any kind means that 
which configuration a pinned manifold stops in depends 
on details of its local dynamics and on its history, in 
a way that cannot occur in the absence of incrtial ef- 
fects. One particularly interesting consequence of this 
is the coexistence in two identical samples at the same 
value of the driving force of a static locally stable config- 
uration, and a moving configuration that will "overtake" 
the static configuration. What the macroscopic conse- 
quences of this are is the primary subject of this paper. 
As stress overshoots can occur more readily than local 
overjumps caused by inertia, we will generally refer to 



both types of effects as stress overshoots, although, with 
slight inconsistency that we trust will not be confusing, 
we will characterize their strength by a parameter that 
we denote M. 

What is the nature of the dcpinning transition as M is 
increased from zero? There are various scenarios one can 
readily envisage. For large enough M, in the model we 
introduce in the next section, an infinitesimal increase in 
F from any pinned state will result in a non-zero aver- 
age velocity at long times. This is because sufficiently 
large stress overshoots always induce other segments to 
move when triggered by an initial segment that moves 
in response to an increase in F. The increased triggering 
will spawn further motion, despite the fact that the stress 
overshoot is only temporary. This process will run away; 
and the manifold will acquire a non-zero average velocity. 
While it might thus seem likely that the transition will 
become "first-order" for large enough M, does it do so 
for arbitrary small Ml If there are regimes in which the 
dcpinning is indeed discontinuous in some way, is there 
macroscopic hysteresis? In what sense? 

More generally, what happens to the depinning transi- 
tion beyond the dissipative limit? If it remains critical — 
at least in some respects — for a range of M, what is its 
nature? Can the quasistatic behavior persist macroscopi- 
cally for small M in spite of the presence of additional mi- 
croscopic hysteresis? If so, what is the size dependence of 
the hysteresis and related phenomena? The system size 
dependence is particularly relevant for geological faults 
for which the statistics of the earthquakes are affected 
both by the nature of the drive and the distribution of 
the "sizes" of faults. 

Several recent papers have undertaken some prelimi- 
nary studies of the effects of stress overshoots in models 
of depinning of elastic manifolds. Reference |lj studied 
several one-dimensional models with long-range elastic- 
ity and stress overshoots motivated by p lanar crack fronts 
driven by applied loads. Reference [Q introduced a par- 
ticularly simple model with short-range static elasticity 
and analyzed its infinite-range limit in which all segments 
of the manifold are equally coupled to each other. In this 
limit, the spatial properties of the manifold are averaged 
away; only time dependence remains and mean field the- 
ory becomes exact. Such mean field models were the 
starting point for theoretical understanding of the finite- 
dimensional physics in the quasistatic limit . Whether 
or not they provide a useful starting point beyond the 
dissipative limit, is one of the questions that we must 
address. 

In this paper, we investigate numerically and phe- 
nomenologically the finite-ran ge version of the stress 
overshoot model introduced in [ft2f. 



2 



A. Outline 

The remainder of this paper is organized as follows: In 
the next section we introduce the basic lattice model on 
which we focus. In Section III the general scaling picture 
is introduced and known results for the dissipative case 
are summarized. In Section IV the critical behavior in 
the moving phase is studied, initially for the dissipative 
case, and then in the presence of stress overshoots. We 
summarize a variety of evidence that the critical behav- 
ior which occurs as the driving force is decreased until 
the system stops is in the same universality class as the 
dissipative case. 

Section V turns to the key aspect of overshoots: hys- 
teresis. We analyze the hysteresis loops that occur 
when the system is stopped from the moving phase and 
restarted by gradually increasing the force. Various puz- 
zling aspects of the data are discussed and some under- 
standing of the hysteresis loops in terms of a low density 
of segments that can be readily triggered by an increase 
in F is reached. In the following section, studies of the 
dynamics and statistics of avalanches that occur as the 
force is gradually increased are presented. These again 
lead to puzzling dependence on overshoot and pinning 
strengths although some aspects of the avalanches ap- 
pear very similar to those in the dissipative limit. The 
data suggest various subtle crossovers may be occurring. 
In Section VII the dynamics of the nucleation of restart- 
ing after a system has been stopped from the moving 
phase are analyzed. It is found that over a substantial 
range of sizes, bubble-like nucleation can occur. 

In Section VIII the puzzling aspects of the various sets 
of data are tentatively resolved in terms of a crossover 
as a function of length scale and system size that mani- 
fests itself in different ways as the various parameters are 
varied. Finally, in Section IX, the conclusions are sum- 
marized and applications to the dynamics of earthquakes 
are discussed briefly, 

In the main body of this paper, we restrict consid- 
eration to weak enough overshoots that they do not to- 
tally change the local dynamics. But for sufficiently large 
M, the overshoots cause dramatic changes in the macro- 
scopic behavior. Although these are interesting, they are 
probably peculiar to certain aspects of the model; results 
on these will be presented elsewhere 13 1. 



II. MODEL 

Near the depinning transition, the dynamics is very 
jerky with segments of the manifold spending most of 
their time stationary or almost so, but occasionally get- 
ting unpinned by the forces from other segments and 
moving forward only to get pinned again by a combi- 
nation of the newly explored random forces and the elas- 
ticity. The inherent discreteness of these local jumps sug- 
gests that we model the manifold as a large number of 



segments that can jump discontinuously from one pinning 
position to another; this is also convenient for numerical 
studies. We define h(x, t) to be the single-valued scalar 
displacement of the manifold from some undeformed ref- 
erence configuration with both the position, x and the 
time, t, taken to be discrete. Note that by constrain- 
ing the displacement field to be single- valued, we exclude 
"overhangs" as well as defects such as dislocations that 
could otherwise occur in periodic systems. The forces on 
a segment of the manifold consist of three terms: the ap- 
plied force F, a static random pinning force rj(x, /i(x, t)), 
and the stress caused by the elasticity <r(x, t). 

The stress depends linearly on the displacements of 
other parts of the manifold via 



<t(x, t) = J2 J2 J *y( T ) h (y> t-r)- Mx, *)> 



y t>0 



where 



(2) 



(3) 



and the sum is over nearest neighbors of x. To model 
stress overshoots, we assume the simplest possible form: 
that the overshoot only applies to neighbors and only 
lasts for one time step so that 



^xy(r 



0) 



1 



(1 + M) 



and 



J, 



M 

xy(l) = - Y > 



(4) 



(5) 



with Z the number of nearest neighbors. With this stress 
transfer, the jump of any nearest neighbor of x induces an 
extra temporary stress on the xth segment. When M — 
0, the stress transferred to the xth segment is simply 
proportional to the static curvature at x. This stress will 
not decrease with time as long as the xth segment does 
not move and the other segments only move forward or 
remain at rest; this limit is thus the dissipative dynamics 
already studied extensively However, for positive M the 
stress on the xth segment caused by a jump forward of 
one of its neighboring segments will first increase by a 
larger amount than in the absence of M and then decrease 
at the next time step to reach its quasistatic value. 

Modeling of the local pinning forces also involves sub- 
stantial arbitrariness. In Ref. we chose randomly 
spaced pinning positions for each segment with uniform 
pinning strengths, but this choice did not affect substan- 
tially the mean field behavior [fl3|| . For our present pur- 
poses, it is more convenient to choose the pinning posi- 
tions for each segment to be uniformly spaced but with 
their yield strengths, the maximum force they can sus- 
tain, randomly distributed. In particular, we take the 
distribution of these yield strengths to be uniformly dis- 
tributed from [0J™ ax ]. 
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Because these random forces pin, or hold back, the 
manifold the corresponding forces ?j(x, h(x)), take on 
negative values in the range ( — / p (x, /i(x)), 0) . With 
this form of the pinning forces, the equation of motion is 
simply given by 

h(x, t+l) = h(x, t) + 9[cr(x, t) + F - / p (x, h(x, t))], 

(6) 

where is the unit step function. The theta function is 
imposed so that a segment can move only forward and 
does so when the net force /(x) on it (the argument of 
the theta function), is positive; otherwise it remains sta- 
tionary. With this dynamics, when a segment jumps its 
displacement always increases by one. As long as v <C 1, 
the (artificial) upper limit on the velocity, this automaton 
dynamics mimics the continuous time motion reasonably 
well. 

Note that in the absence of elasticity, with M = 0, 
F c = f™ ax with each segment becoming stuck on an 
anomalously strong pinning site. In the presence of the 
elasticity, not all of the segments can be simultaneously 
pinned on strong pinning sites, and the critical force will 
decrease. For weak pinning forces, however, the defini- 
tion of the applied force in this model is somewhat patho- 
logical: when released from a pinning segment, a segment 
can jump forward far enough that the total force on it 
becomes negative and it can then pull forward other seg- 
ments resulting in overall motion even if F is negative. 
To make it more realistic, one could replace F by F + 1; 
so that there is always enough force to make the forces 
at pinning segments non-negative. With non-zero M, an- 
other adjustment should really be made as more realistic 
forms of stress overshoot involve a concomitant negative 
force on the segment that has moved. Since the zero 
of the applied force is entirely a convention, we will not 
make these adjustments; this will mean, however, that in 
some regimes the critical "force" will be negative. 

Several additional aspects of the model need to be 
specified: the initial conditions and the order of the up- 
dating. To avoid lock-step or other "faceting" like behav- 
ior [Q , we choose the pinning positions of each segment 
to be offset from one another by random amounts in the 
interval [0,1]. Because of the integer character of the 
jumps, the fractional part of the displacement of a given 
segment does not evolve with time. It might be thought 
that the displacement-independent randomness induced 
by this constraint could dominate over the randomness 
of interest, especially as far as determining the variations 
of critical forces, etc. in finite size systems. [Indeed, just 
such an effect does occur for systems with periodic ran- 
domness such as charge density waves ||. But in the 
present case it can be shown that the additional ran- 
domness is analogous to a spatially random force that 
is the derivative of a random function. This, combined 
with the statistical tilt symmetry, mean that its effects 
are subdominant for large systems (although they could 
give rise to additional corrections to scaling.) The fact 



that the variations in the critical forces in finite size sys- 
tems decrease substantially faster than the inverse square 
root of their area supports this assertion. 

The updating of the displacements are done in paral- 
lel after computing all of the stresses. While there are 
alternate sequential methods of updating, this parallel 
method requires the least amount of computation and 
does not appear to introduce any troublesome artifacts 
in the regime of smallish M of primary interest here. 

Finally, to limit boundary effects, we impose periodic 
boundary conditions on h(x, t). This is especially impor- 
tant as we will use finite size scaling to analyze much 
of the data; this is substantially more straightforward 
with periodic boundary conditions. Our simulations are 
restricted to two dimensions which we chose because of 
the availability of the widest range of system sizes with- 
out running into the complications associated with very 
large stresses that arise in one dimensional depinning 
with short range interactions. We study systems of size 
L x L up to 256 x 256 with most of our "large" system 
data on 128 x 128 samples. 

III. SCALING AND DISSIPATIVE DYNAMICS 

Before presenting new results for the systems of in- 
terest with stress overshoots, we briefly summarize the 
scaling behavior that obtains near the critical force in 
the absence of stress overshoots; i.e. for M = 0. As the 
force is adiabatically increased from zero, local instabili- 
ties lead to a succession of avalanches, most of which will 
be small, but which can occasionally become large as the 
unique critical force F c q is approached. Above F c q the 
mean velocity in the statistical steady state rises contin- 
uously with an exponent (3. The motion is jerky out to 
length scales of order the velocity correlation length 
which diverges at the critical force as 

~ (F-F c0 y ■ (7) 

The characteristic time for relaxation on scales of order 
6« is 

r~e v (8) 

and in this time the manifold typically moves forward by 
an amount 

Ah ~ e (9) 

These three exponents, v, z and £ characterize the scaling 
behavior near the transition. The velocity exponent is re- 
lated to these via the observation that the mean velocity 
is of order the characteristic displacement per character- 
istic time so that 

V ~ ~ (F - (10) 

with 
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= {z-Qv. 



(11) 



In the pinned phase, the critical behavior as the force 
is adiabatically increased can, in the absence of stress 
overshoots, be related to that in the moving phase re- 
viewed above. In particular, the scaling of the dynamics 
and shape of the avalanches, the probability that they 
will be large, the divergence as F c q is approached of the 
cutoff size in their distribution, and the "roughness" of 
the manifold at the critical point are all given in terms 
of the same three exponents. In Section VI we will dis- 
cuss the avalanches in detail, but for now we focus on the 
macroscopic behavior such as the velocity in the moving 
phase and the mean displacements in the pinned phase. 

The mean displacement in response to a spatially vary- 
ing applied force yields, via a statistical symmetry of the 
system, a scaling law that relates two of the exponents. 
This relation can be derived from the average static po- 
larizability 



These yield 



(17) 



(18) 



Recently Chauve et al [y_5[ have computed these expo- 
nents to second-order in e, obtaining 



and 



C ~ -(1 + 0.14331e) 



z re 2- -e- 0.04321e^ , 



(19) 



(20) 



although there are some doubts about the validity of 
these second-order results M. 



X(q,w = 0) = 



dh{q,LU = 0) 



de(q) 



(12) 



to a perturbing force e(q) cos(q ■ x). A change of variables 
to h(x,t) — h'(x,t) — V~ 2 e(x), yields an equation of 
motion for h! that is statistically identical to the original 
one for h, independent of the perturbing force. Therefore 



dh(q,w = 0) 



de(q) 



l-J(q) <? 2 ' 



(13) 



Since the polarizability should scale as Ah/AF ~ 
' v X(q£), with X a scaling function, this yields the 
scaling law 



C + ± = 2. 

V 



(14) 



A. Dissipative exponents 

The critical exponents for M = take simple mean 
field values oi i> — 1/2, z = 2, characteristic of diffusive 
dynamics, and C = above the critical dimension of d c = 
4 for short range elasticity. The velocity scales with (3=1 

!• 

Below four dimensions, renormalization group expan- 
sions have been performed that justify the scaling laws 
and claims of universality as well as yielding results for 
the exponents as expansions in powers of 



= 4- d 



(15) 



(16) 



and 



IV. CRITICAL BEHAVIOR IN MOVING PHASE 

We will shortly turn to presentation of our numerical 
results for the critical behavior in the moving phase. But 
first, it is instructive to summarize the behavior found in 
the mean field limit for small M and to consider sev- 
eral possible scenarios that might obtain in short-range 
systems. We can then determine which scenario is most 
consistent with the data. 



A. Mean field limit and scaling scenarios 

In the mean field limit all of the sites are coupled to all 
of the others so that the number of "nearest neighbors" 
Z is equal to the number of segments, TV (more precisely, 
Z = N — 1). In this limit, the critical force is found to 
be unchanged for M less than a critical value, M c = 1. 
The velocity versus force curve is modified for any non- 
zero M, however, but the exponent (3 remains at its qua- 
sistatic value of j3 = 1 for M < M c . The other universal 
properties of the transition are also unchanged for small 
M, including the lack of hysteresis in steady state and 
the asymptotics of the distribution of large avalanches as 
the critical force is approached from below. 

The simplest scenario for the short range systems of 
interest would be like that of the mean field limit: un- 
changed critical behavior and no macroscopic hysteresis 
for small M. But previous work has shown that this 
cannot be the case: As shown in references |jTlf , |]l6| , any 
stress overshoot will cause the critical force to be shifted 
downwards, in the sense that dynamic behavior can per- 
sist for some (M -dependent) range of forces below the 
quasistatic critical force F c q. This implies that some form 
of macroscopic hysteresis can exist since locally stable — 
at least linearly stable — static configurations exist up 
to F c q. But whether such configurations are non- linearly 
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stable to, for example, an arbitrarily small increase in F, 
is a question of substantial importance to which we will 
return later. For now, we focus on the moving states and 
how they stop as the force is lowered. 

The simplest scenario that cannot immediately be 
ruled out is a modified version of the mean field scenario: 
a velocity versus force curve with a well defined critical 
force, F}, that is non-hysteretic as the force is decreased; 
history independent steady states above F*\ and critical 
exponents, (3, , unchanged from their quasistatic val- 
ues. In renormalization group language, this would cor- 
respond to M being an irrelevant perturbation at least 
as far as behavior in the moving phase. 

If this mean-field-like scenario indeed applies for small 
M to the finite-range model, we expect the following 
scaling behavior for both M and the proximity to the 
M-dependent critical force, 

f = F-F c l (M). (21) 

small: 

v(M,F)~fB(^) (22) 

with the crossover exponent 4> < indicating the irrele- 
vance of M, and B a scaling function. In analogous situa- 
tions in equilibrium statistical mechanics, if a parameter 
such as M does not change the nature of the transition, 
the effect of M on F* can be taken into account per- 
turbatively. Because of the singular nature of the critical 
fixed point that describes the quasistatic depinning — re- 
sulting, in part, from the absence of thermal fluctuations 
but, more essentially, from the jerky nature of the motion 
— the "analytic parts" might not be smooth functions of 
M, but under the assumption that they are, we expect 
that 

F r l (M) = F c0 - aM - bM 2 + ... . (23) 

Nevertheless, there can be singular corrections to scaling 
determined by the form of B and the crossover exponent 

This scenario, in which M is irrelevant when it is small, 
we call the dissipative scenario. We note, however, that 
this scenario is compatible with a change in behavior at 
a critical value of M as occurs in mean field theory. This 
would give rise, for M close to its critical value, to a 
crossover to some kind of multicritical behavior emerging 
at larger velocities. 

More interesting behavior would occur if the qualita- 
tive mean field results on the effects of small M do not 
simply carry over to the finite-range case. This would be 
the case if M is a relevant perturbation and would cor- 
respond to a crossover scaling function, such as B in Eq. 
(22) with <fi positive. A relevant perturbation M would 
yield a singular correction to the critical force of the form 

[FcO-Fi^Ug-M 1 ^, (24) 



which would dominate over the leading (or subdominant) 
analytic shift if <f> > 1 (or <j> > |). Earlier numerical 
results by Ramanathan and Fisher suggested that this 
might be the case |0|. 

There are two simple scenarios for the velocity versus 
force if M is relevant. One possibility is that the depin- 
ning transition is driven discontinuous immediately for 
any non-zero M. The average velocity would then have 
a discontinuity of 

v mm ~ Ml (25) 

We call this the first-order scenario. 

If M is relevant but the transition is still continu- 
ous, one would expect it to be in a different critical 
universality class. In this case, the scaling function 
B(y — > oo) ~ y p so that v(f) remains continuous but 
with a new exponent 

P = f3 -<j> P (26) 

with /3q the quasistatic value. Asymptotically close to the 
critical force the new critical behavior would obtain but 
for / > > M * , the average velocity curve would crossover 
to the dissipative behavior. We refer to this as the new 
universality class scenario. 

While more exotic scenarios may be possible, we will 
limit our consideration to the three scenarios enumerated 
above. 

Note that we have explicitly not considered scenarios in 
which the velocity in the moving phase is hysteretic. Al- 
though we cannot rule this out entirely, the fact that the 
random environment through which the manifold moves 
acts, to some extent, like thermal noise, suggests that 
if there were more than one possible moving phase for 
the same applied force, there would be some stochastic 
process by which the system could jump from one to the 
other. The result of this would be that, as in equilibrium 
transitions, true "coexistence" would not be possible over 
a range of parameters. This argument does not apply to 
coexistence between static and moving phases, as the for- 
mer are not subject to time dependent "noise" . 

B. Numerical results in moving phase 

Our primary numerical results in the moving phase 
were carried out on square two-dimensional samples of 
linear dimension L = 128. The maximum pinning force 
/™ aa: , was chosen to be either 0.5 or 1.0 so that the typical 
pinning force is comparable to the change in elastic forces 
caused one neighbor of a segment jumping. 

We focus on the results for M small enough that sub- 
lattice effects are not too important — for /™ aa: = 0.5, 
we study M < 0.8. Fig. 0(a) shows v(F) for < M < 0.6 
in increments of 0.2. To generate these curves we start, 
for each M at F > F*, with almost flat initial condi- 
tions — segment displacements random in (0, 1) — and 
then decrease the applied force very slowly until Fj-(M) 
is reached. 
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FIG. 1. (a) Steady state velocity, v(F), for L = 128 and 
= 0.5. On the actual plot, as well as on all remaining 
ones, f™ ax is denoted by fpm. The applied force is changed by 
an increment of 0.005 after the manifold has been equilibrated 
for 12,000 time steps (see text). Initially, F > F}(M) and 
then it is decreased until v = before being increased back 
up again to its starting value. The filled arrows indicate this 
cycling of F. The errors bars are smaller than the symbols, (b) 
The same as (a) but with the curves shifted by an approximate 
Fi(M). 



For our finite systems with L = 128, Fj-(M) is defined 
as being the force below which the system halts after 
12, 000 time steps, a value chosen so as to be long enough 
for transients to decay, but not so long that rare con- 
figurations of the randomness with anomalously strong 
pinning forces can dominate. 

For M — 0, there is a finite-size crossover regime 
in which the system may stop due to an anomalously 



strongly pinned region and the infinite system behavior 
will no longer be observed. This tends to occur when 
the average velocity while it is still moving is of order 
L~~ . For L — 128, this value is about 0.013 and with 
this average velocity, the manifold will typically travel 
a distance several times the characteristic displacement 
i> within 12,000 time steps. As we shall see, our esti- 
mates are that at least one measure of the characteristic 
time scale, which grows as L z , is only weakly dependent 
on M . Thus, we use this same criterion for non-zero M 
while being aware that it may bias our scaling results 
slightly such that if time scales do change substantially 
with M we may be observing more of either finite-size 
effects or non-equilibrium effects, the former if the time 
scale decreases with M and the latter if it increases with 
M. For other system sizes, the equilibration time is de- 
creased correspondingly, roughly with L 2 °, i.e., according 
to the dynamic scaling found in the dissipative limit. 

We first study the dependence on M of the critical 
force, F^, below which the steady state motion ceases, in 
particular to test whether Fj_(M) is a singular or smooth 
function of M. See Fig. g The results, along with a 
quadratic fit, 



F c x (0) - Fj(M) = 0.27(±0.01)M - 



0.012(±0.01)M 2 
(27) 



are shown in Fig. |[ [Note that if we had included a 
constant in the fit, the constant would have vanished 
within one standard deviation, as it should, and the lin- 
ear coefficient would have been only slightly modified to 
0.26 ± 0.01.] A natural expectation — although overly 
naive — is a linear decrease of F} by an amount M/Z. 
For small M, this appears to work rather well. The rea- 
son for this linear shift and the corrections to it will be 
discussed later. 

The analytic fit should be compared to a fit — with the 
same number of parameters — to an arbitrary power-law 
M dependence of the shift in Fj such as would obtain if 
M were a relevant perturbation. The best fit to the F* 
data yields an exponent 4 = 1.18 ± 0.02; note, however, 
that by eye the quadratic fit looks slightly better than the 
power-law fit. Although with a weakly relevant M with 
a crossover exponent less than unity, as the power law 
fit suggests, one would presumably have a linear analytic 
term as well and thus the inferred cf> w 0.85 should not 
be taken too seriously, In any case, one must ask whether 
it is consistent with our other data. It does not appear 
to be: if we use this value of 4> to try and find a scaling 
function B(jj) for the velocity data of Fig. [3](a), the 
curves do not collapse. This suggests that either M is 
irrelevant, or that it is sufficiently weakly relevant that 
the crossover exponent <f> is small enough that it would 
not dominate the shift in . Other data, as summarized 
below, suggests that, in fact, M is irrelevant, at least for 
the steady state moving phase. 
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FIG. 2. Critical force, F r l (M = 0) - F}(M), as a func- 
tion of M indicated by filled circles; the solid line is the re- 
sult of a quadratic fit bM + cM 2 with b = 0.27 ± 0.01, and 
c = 0.12 ± 0.01 while the dotted line is the result of a fit to 
aM 1+D with a = 0.384 ± 0.003 and D = 0.18 ± 0.02. The 
error bars are smaller than the symbols. 



The mean velocity data can be used to obtain the criti- 
cal exponent (3 and see whether it depends on M . Figure 
shows a log- log plot of the v(F) curves with the best 
F}(M) value, that which makes the curve the most lin- 
ear, determined by hand for each M. The error bars 
indicated are the rms variations in the average velocity 
over 10 samples. The values of the applied forces used 
are separated by an interval of 5 • 10 -4 , one-tenth of those 
used in Fig. 0(a). 
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FIG. 3. Log-log plot of the average steady state v as a 
function of / = F — F^{M), where F}(M) is varied to op- 
timize the linearity of each curve within the scaling regime. 
The symbols used are the same as in Fig. |l|(a) , with the filled 
squares representing M = 0.2, f™ ax = 1.0. Each curve repre- 
sents an average over 10 samples with the error bars represent- 
ing the rms sample-to-sample variations in the velocity. The 
applied force increment is 0.0005 for approximately v < 0.1, 
while for v > 0.1 the force increment is 0.005. For the weaker 
randomness, the best fit slope for M = 0.0 is 0.66 ± 0.03; for 
M = 0.2 it is 0.66 ±0.02; for M = 0.4 it is 0.68 ±0.05; and for 
M — 0.6 the slope is 0.67±0.10. For the stronger randomness 
at M = 0.2, (3 = 0.66 ±0.02. 



The velocity critical exponent (3 inferred from these 
data is, for M = 0, 



f3 = 0.66 ±0.03. 



(28) 



Surprisingly, this value appears to be consistent within 
one standard deviation with the data for all the M values 
shown. Even for M — 0.6, we find 



(3 = 0.67 ±0.1 



(29) 



although the straight line fit is only over one and a half 
decades in the reduced force, /, substantially less than 
the three decades of the fit for M — 0. From Fig. |l|(a), it 
is apparent that this reduced range of scaling is primarily 
due to a larger amplitude for the singular velocity for 
larger M. See Fig. 0(b). 

The data suggest that the evidence is at least con- 
sistent with the dissipative critical behavior obtaining 
asymptotically for small M; i.e., with M being an ir- 
relevant perturbation. 

We must be careful, however, especially as Figures 0(a) 
and (b) , and indicate that the minimum average veloc- 
ity is increasing with increasing M . Might this obser- 
vation suggest that there is a discontinuity opening up 
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as a power of M, perhaps suggesting that the transition 
is driven discontinuous immediately? We can check the 
data against the finite-size crossover behavior expected 
in the M — quasistatic limit; this would yield 

v min ~L~~. (30) 

Figure f| tests for this scaling and finds it to be consistent 
with the data: even for M — 0.4, the slope of the log-log 
plot is — 0.91±0.01, which agrees within one standard de- 
viation with the M = 0.0 slope. While we do observe an 
increase in the minimum average velocity with increas- 
ing M at fixed size, it is not the size-independent power- 
law increase with M that would have been expected if 
the transition became discontinuous. Instead, there is an 
M-dependent coefficient 
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FIG. 4. Plot of the minimum spatially averaged velocity 
Vmin as a function of both M, L and f™ ax . For the determi- 
nation of v m in, see the text. The symbols for the various M's 
are the same as in Fig. |l](a), except for the left triangles which 
represent M = 0.1. The best fit slope for the M = 0.0 curve 
is -0.90 ± 0.02, for M = 0.1 it is -0.89 ± 0.03, for M = 0.4 
the slope is -0.91 ± 0.01, and for M = 0.6 it is -0.83 ± 0.04. 
Each set of data are averaged over 10 samples. 



associated with the average velocity, just as occurs for 
small M in mean field theory. Similar minimum veloc- 
ity data was also obtained for stronger randomness, with 
jrnax _ ^ ^ nug gee faa,t the finite-size data are con- 
sistent with the dissipative scenario, with small stress 
overshoots being irrelevant for the velocity versus force 
curves. 

But it is still possible that a new critical universality 
class is emerging for small M if p << 1, so that the 
emergence of the new universality class would be diffi- 
cult to detect by simply measuring the velocity exponent 



as this would be little changed from its dissipative-limit 
value. We therefore look more closely at the finite-size 
crossover regime to investigate whether other aspects of 
the behavior really look similar to the M = quasistatic 
dcpinning. 

Anticipating that it might be the stopping behavior 
that would distinguish between quasistatic and overshoot 
dynamics, we have explored some of the dynamics of the 
stopping process. Given that there is a distribution of 
Fj's, we let the manifold equilibrate at an F three stan- 
dard deviations above the average of Fj(M). We then 
lowered F to the average of F*(M) and waited for the 
manifold to come to a stop. Figure || plots the instanta- 
neous spatially averaged velocity v(t) averaged over many 
samples as they come to a stop at time t stop . Here we ex- 
plicitly see that the various M samples come to a stop in 
a similar gradual manner; only the amplitudes vary. The 
differences between the curves are presumably due to the 
M dependence of the amplitude of the steady state v(F) 
curves. The final stages of the stopping process can be 
analyzed — at least for the dissipative (M = 0) case — 
by a simple scaling argument. Once only a small fraction 
of the system is still moving, the average velocity will be 
inversely proportional to the area L d . Since the velocity 
scales as and lengths as t* , we expect 



v(t) 



stop 



t) 



(d+(-z)/z 



L d 



(31) 



the data on a log-log plot, Fig. ||, are reasonably consis- 
tent with this for all values of M tested. 
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FIG. 5. Log-log plot of v(t Btop — i) as a function of the 
time before the motion ceases at i s to P - The data were ob- 
tained from a uniformly random initial condition on the inter- 
val [0, 1] and equilibrated at an applied force higher than F} 
by three times the standard deviation of this critical force. Af- 
ter equilibration, the applied force was lowered to the average 
F}, whereupon v(t) is recorded until the manifold eventually 
stops. The solid line is the theoretical expectation within the 
scaling regime. The slopes of the numerical data, in order of 
increasing M, are 0.831 ± 0.001, 0.846 ± 0.002, 0.981 ± 0.003, 
1.085 ± 0.004, and 0.975 ± 0.004 respectively. Note that for 
the M = 0.8 curve, f™ ax = 1.0. Each curve represents an 
average over those of approximately 500 samples. 



FIG. 6. The average maximum width W m „ of stopped 
manifolds as a function of L for various M's and /™ a:E 's. The 
data are taken after the manifold has stopped moving follow- 
ing a gradual decrease of the driving force. Each W ma x is an 
average over 10 samples. The slopes can be used to extract 
the exponent For the left part of the figure, f™ ax = 0.5; 
the M = 0.0 data are best fit by 0.75 ±0.05; M = 0.2 yields a 
slope of 0.72 ± 0.03; and for M = 0.6 the slope is 0.75 ± 0.05. 
For the the right plot, f™ ax = 1.0; for the M = 0.0 curve, the 
slope is 0.75 ±0.02; for M = 0.2 it is 0.81 ±0.02; the M = 0.6 
slope is 0.75±0.01; while for M = 0.8, the slope is 0.79±0.01. 
Note that the overall magnitude of Wmax is slightly larger for 
the stronger pinning. 



We can also probe the finite-size crossover regime in 
terms of the roughness. We define the maximum width of 
the manifold IFmax as the absolute value of the maximum 
deviation of the displacement from its spatially averaged 
value, h: W max =< max|/i(x) — h\ >. This should scale 
as with £ « 2/3 for M = 0. If M is irrelevant, then 
the same should hold true for all small M. Figure ^ 
demonstrates that the maximum width does indeed obey 
this scaling with system size, but with apparent values 
of Q that are somewhat larger than | even for M = 0.0 
for which ( = 0.75 ± 0.05 is inferred from the data with 
jmax _ q 5 Although the exponent appears to be M- 
independent, there is an M-dependent coefficient with 
the overall width increasing — albeit only slightly - 
with increasing M . For the stronger randomness data, 
jmax _ ^ this tendency is less strong but that is because 
the value of M at which the behavior changes character 
is larger for stronger pinning. Overall it appears that, at 
lpa&t up until M = 0.6, the finite-size crossover regime 



Iodks similar to the dissipa:;ive M = case 
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We can also determine £ from the power spectrum of 
the displacement just after the motion has stopped from 
the moving phase; we use the same equilibration time 
and applied force increments as in Figs, ^(a) and [7] (b). 
In two dimensions, 



< \h(k)\ 2 >= / d 2 x e* k x < h(x)h(0) > ~ k 



-2C-2 



(32) 



where the brackets denote averaging over samples. From 
Fig. ||(a) we see that fitting to this form over a range of 
one and a half decades for samples with L = 256 yields 
C = 0.72 ±0.02 and 0.74 ±0.01 for M = 0.0 and M = 0.5 
respectively, again in mutual agreement and similar to 
the values from the maximum width discussed above. 
But again, apparently slightly larger than the C = § value 
from the first-order epsilon expansion. Again, the cumu- 
lative evidence suggests that at least up until M > 0.7, 
it is not likely that a new universality class is emerging 
in the "equilibrium" moving or stopped phase as v — > 0. 

It is useful to compare our values of the exponent £ 
with those previously obtained. Leschhorn et. al. found, 
in numerical simulations, ( — 0.75 ± 0.02 for M = 0.0. 
Jl7t The most solid theoretical result is that | is a lower 
bound for (. This comes from application of finite size 
scaling to the connections between the variations of the 
critical force and the correlation length exponent — see 
discussion below. While Narayan and Fisher had ar- 
gued that the value of | is exact to all orders in e, Chauve 
et. al. jl5| have computed the exponents to second-order 
in e and found that the roughness exponent is increased 
to C = f (1 ± 0.14331e) ± C(e 3 ), which naively extrap- 
olated to two dimensions yields £ « 0.86 substantially 
higher than the value inferred numerically. Whether or 
not this discrepancy is due to the neglect of terms higher 
order in e, to some problem with the expansion, or to 
corrections to scaling remains in doubt. It is worth not- 
ing in this context that for one-dimensional systems with 
the long range interactions appropriate for crack fronts, 
Ramanathan and Fisher [jll| found that to obtain reli- 
able and universal values of £ analysis of corrections to 
scaling were needed. With these included, they found a 



10 



value of £ very close to that predicted from the first-order 
e-expansion without any higher order corrections. 
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FIG. 7. (a) Log-log plot of the square root of the power 
spectrum of a manifold in steady state both just above the 
depinning transition, open symbols, and just after the mani- 
fold has stopped after the force is decreased gradually, filled 
symbols. The various sets of data have been shifted along 
the vertical axis for clarity. For the moving configurations, 
F = -0.0775 for M = 0.0, F = -0.1389 for M = 0.2, and 
F — —0.2436 for M = 0.5. For each M the data are averaged 
over 100 samples of size of L = 256. The roughness exponent 
£ can be extracted from the linear fits over the appropriate 
wave vector interval: the expected slope is — (1 + f). The 
quoted error bars are a measure of the dependence of the ap- 
parent slope on the range over which the data are fit. The 
slopes of the moving data are —1.65 ± 0.01, —1.65 ± 0.02, 
-1.69 ± 0.01 for M = 0.0, M = 0.2, and M = 0.5 respec- 
tively. While the slopes for the stopped data are —1.72 ±0.02 
and -1.74 ± 0.01 for M = 0.0 and M = 0.5. Note that at 
the smallest wavevectors flattening is evident in the moving 
data but not the stopped data, as expected. Also note the 
appearance of a peak at the zone corner in the moving data; 
this represents the tendency for alternating sublattice motion 
that occurs, especially for larger M. (b) We plot data for 
M = 0.8, fpm = 1.0, and L = 256 for both stopped (light 
dots) and moving (dark dots) configurations We extract a C, 
of 0.75 ± 0.01 for both curves. These curves have not been 
shifted with respect to each other. We note that the data 
in both (a) and (b) has been thinned out for clarity, (c) 
Same plot as (a) without any shifting of the data and over 
a narrower region in k-space so that the amplitudes can be 
compared. 



For completeness, we have also studied the spatial 
power spectrum of the displacements just above the de- 
pinning transition. Figure ?](a) and (b) shows the square 
root of the averaged power spectrum for v « 1. The 
roughness exponents obtained for M = 0.0, 0.2, 0.5 are 
again roughly in agreement with each other — £ = 
0.65±0.01, 0.65±0.02, and 0.69±0.01 respectively — but 
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slightly smaller than those found for the stopped man- 
ifold. In the moving phase, however, there should be a 
crossover at long wavelengths, k£ v ~ 1, to the Edwards- 
Wilkinson universality class with only logarithmic rough- 
ness |Q . This arises because on scales larger than , the 
motion makes the randomness appear like white noise in 
both space and time and the displacement correlation 
function becomes 



for all sufficiently small M as the transition is approached 
from above; this despite the A/-dependent shift in F£ and 
the concomitant hysteresis that is possible because of the 
existence of linearly stable static configurations up^lp F c q 
which is I larger than F* . — We will analyze this para dox 
later. For now^jfe appears that the mean-ficld-like sce- 
narifcOnafs- won out over the first-order scenMiojand the 



new universality class scenario 



< h(k,uj)h(-k,-Lu) >r 



1 



Dk 2 



(33) 



with D(F) an effective diffusion constant. This crossover 
is observed for the smallest k in Fig. (?] (a) where the slope 
of the power spectrum decreases in contrast to the data 
taken after the motion has stopped. We note that this 
difference in slope between the moving and stopped spec- 
tra is not as prominent when the randomness is stronger 
as shown in Fig. (b). 

For the moving configurations we observe a peak in the 
power spectrum at k = (71", w). This peak is caused by the 
tendency of one segment's motion to trigger jumps of its 
neighbors at the next time step. The structure and am- 
plitude of the peak looks similar for all of the M's shown, 
although the wavevector dependence indicates that there 
are somewhat more segments participating in the sublat- 
tice behavior at M = 0.5 than at M = 0; but only about 
a third more. 

The dynamic exponent in the moving phase, z, can be 
determined in various ways from data near to the critical 
force. We first study the non- equilibrium roughening of P 
the manifold starting with almost flat initial conditions. 
We define 



w 2 (t) = ((h( X ,t)-h(t)) 



2\ _ i-r 



2C 
t~. 



(34) 



with the overbar denoting spatial averaging over the sam- 
ple. The scaling behavior is expected for times short com- 
pared to the critical correlation time r which diverges at 
F* . Once the roughness exponent has been calculated in- 
dependently from the same set of configurations in Fig. 
0(a), z can be extracted from the log- log plot of w 2 (t) vs. 
t as is done in Fig. |(a)(b) For M = 0.0 and M = 0.2 
we see that the dynamic exponents z are very similar: 
z = 1.51 ±0.03 for both. However, for M = 0.5, the data 
are somewhat further above the transition as they cor- 
respond to an equilibrium average velocity of v = 0.08. 
While these data could be fit with the same z over a 
limited range of times, there is clearly some new physics 
emerging: an upward curvature on the log-log plot and 
a substantial — more than a factor of two — overshoot 
in the velocity before it settles down to its steady state 
value. This effect is related to the change in the dynam- 
ical onset of the motion to which we will turn in Section 
VII. 

Aside from the transient effects associated with ap- 
proach to steady state, all our measurements in the mov- 
ing phase suggest that the critical behavior is most con- 
sistent with the dissipative universality class obtaining 
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FIG. 8. (a)Dynamical roughening as the manifold starts 
moving just above from almost flat initial conditions. The 
mean-square width of the manifold as a function of time is 
shown on a log-log plot. The curves has been shifted along 
the vertical axis for clarity. The expected slope in the scal- 
ing regime of these curves is — so that we can extract the 
dynamic exponent z after measuring £ independently (Fig. 
0(a)). The linear fits yields a slope of 0.86 ±0.01 for M = 0.0 
and 0.86±0.01 for M = 0.2. The scaling regime should obtain 
only for t << £ z ; since we are above the transition £ << L so 
that a crossover to time-independent roughness should occur 
when t ~ £ z . The crossover scale is in rough agreement with 
theoretical estimates. For M = 0.5 the behavior is rather 
different that smaller M: they exhibit too much curvature to 
be readily explained by corrections to scaling and a substan- 
tial overshoot which we believe is due to the merging of many 
nucleating bubbles as the steady state is approached. The ap- 
parent slope in this regime is 1.20±0.02, although this is most 
likely just an effective exponent. For the stronger randomness, 
f™ ax = 1.0, at M = 0.8 a fit would yield 3£ = 0.99 ± 0.01. 
Despite this larger dynamic roughening exponent, the static 
roughness exponent is also larger and the inferred z is not 
much changed from the M = data, (b) Same as (a) but 
without any shifting of the curves. 



C. Amplitude ratios 

If we accept that M is irrelevant for the steady state 
moving critical behavior, the observed changes as M is 
increased in the moving phase for F just above Fj appear 
to be primarily attributable to the increased amplitude 
A V (M) of the velocity near F c x , v ~ A v (M)f f3 . As seen 
in Figs. 0(b) and ||: For larger M, the velocity rises more 
rapidly with increasing F until it becomes close to | at 
which point sublattice effects set in. A similar increase 
in the amplitude of the velocity is observed in mean field 
theory, with the exponent (3 unchanged for small M but 
the amplitude of the velocity growing with M. In both 
cases, this implies that the width in F of the depinning 
transition narrows with increasing M — naively just by 
the narrowing of the range of F — F^ over which v is 
small. 

If the universality class of the critical behavior in the 
moving phase is independent of M over a range of M, 
then there should be various universal relations between 
non-universal coefficients: universal amplitude ratios. A 
priori, we would expect three non-universal scale factors 
associated with the scaling relationships between length 
and, respectively, deviation from criticality, F — Fj; dis- 
placement, and time. We can define these, for example, 
by the scaling of the correlation length, 

6.»(^)"i <»> 

mean square displacements at separations smaller than 



([h(x,t)-h(0,t)f)KA 2 h \x\ 2 <; (36) 
and velocity 

A t * v 
. Ah( F-F} \ p 
~ A\ A F J 

= A V {F-F^f (37) 

where we measure all lengths in units of the lattice con- 
stant. But in the absence of stress overshoots the statis- 
tical "tilt" symmetry of the system that relates the expo- 
nents £ and v via the triviality of the averaged response 
to a static spatially varying additional applied force, also 
relates the associated coefficients: 

^f=C K K (38) 
Ah 

where Ck is a universal dimensionless coefficient and 
K is the long- wavelength elastic constant, which, in our 
model, is simply the inverse of the coordination number 
Z. 

For M = 0, we thus expect that amplitudes of scal- 
ing laws that only involve the three exponents z, v and £ 
should be expressible in terms of the two amplitudes Ah 
and At only. For example, the rms variations in the crit- 
ical force F^ , in finite size systems should be expressible 
in terms of Ap or, via Eq. (36), Ah'- 

^v&r{FHL))^C A C F L-i (39) 

with Ca universal. 

It is not clear, a priori, whether the statistical tilt sym- 
metry argument can be applied in the presence of stress 
overshoots and the concomitant local hysteresis. This is 
because, in essence, it relies on the history independence 
of linear response of at least some quantities. In the mov- 
ing phase, on which we are currently focusing, it seems 
reasonable that the argument should apply and the am- 
plitude ratios hence be related as in the dissipative case. 
But we should remain alert to the possibility that appar- 
ent failure of expected scaling laws may be due to this 
assumption. 

Qualitative examination of the numerical data for the 
roughness (e.g. Fig. 0(a) and the variations in the crit- 
ical force Fj, (see Fig. 0)) suggest that neither of the 
amplitudes Ch and Cf are strongly dependent on M in 
the range studied. One would then guess that the depen- 
dence of the amplitude of the velocity, A V (M), on M is 
primarily caused by a decrease in the characteristic time 
scale as M increases. This will have consequences for the 
behavior of other quantities as will now be discussed. 
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FIG. 9. Log- log plot of the rms variations of F}(M) as 
a function of the system length L. The inverse correlation 
exponents 1/u can be inferred from the slopes. For the 
M = 0.0 curve the slope is -1.29 ± 0.03; for M = 0.2, it 
is -1.29 ± 0.02; for M = 0.4, the slope is -1.18 ± 0.06; and 
for M = 0.6, it is -1.22 ± 0.03. For the stronger random- 
ness data, when M = 0.8, the slope is 1.21 ± 0.03, while for 
M = 1.0, it is 1.19 ±0.05. The dotted line along the M = 0.0 
curve is the result of a fit that includes a correction to scal- 
ing with a correction exponent | and correction amplitude 
B v — 0.78 ± 0.01; the inferred modified correlation length 
exponent is —1.30 ± 0.03. 



A useful quantity to study is the temporal fluctua- 
tions in the instantaneous spatially averaged velocity v(t) 
in finite-size samples: this has information about both 
length and time scales. For a fixed average velocity v, 
the magnitude of the fluctuations should increase with 
increasing M because the system is effectively closer to 
the transition. Since regions of size of order the velocity 
correlation length will fluctuate roughly independently, 
the variance of the instantaneous velocity should, in a 
system much larger than £, be 



Since we have not yet defined the correlation length pre- 
cisely, this could well serve as its definition, thereby fix- 
ing the definition of the amplitude Ap. The universality 
of the amplitude ratios can then be checked by compar- 
ing the M dependence of the amplitude of var[v(i)] with 
those of the variations in F* and the roughness. These 
results are presented in Table I. 
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TABLE I. Universal amplitude ratio data. The parenthe- 
ses denote the uncertainty in the last digit of the quoted 
amplitudes. The first row is M; all the data shown have 
/" OI = 0.5. The second row is obtained from the amplitude 
of the Fourier transform of the h — h correlation function for 
L = 128 slightly above F^ like the data shown in Fig. |^(c) 
for L = 256. The third row is the combination of amplitudes 

Ah/At determined from the combination v( ' ' " — )^, 

where F(0), v, (3 and v are all measured. The fourth row is ob- 
tained from the ratio of the second to the third rows The fifth 
row is obtained from the combination —, -. The 

A h {F-Fl(M))P 

sixth row is the expected universal ratio Ck determined from 
the second, the fourth, and the fifth rows: note the agreement 
within error bars of all the columns. The seventh row is de- 
termined from Af = CkKAh, with the M — 0.0 value for 
Ck used. The eighth row, the velocity amplitude A v , is ob- 
tained from the relation 4^ -nr ■ The ninth row shows the A v 

obtained from the best power law fits of the velocity data of 
Fig. ^. Note the agreement within the errors with the eighth 
row, with the exception of the largest M for which there is a 
small apparent inconsistency. The tenth row is obtained, via 
Eq. (^) from the integral of the ratio of the time-dependent 
velocity- velocity correlations J F(t)dt to r(0) for particular 
choices of the forces close to the critical force. The eleventh 
row, the velocity correlation time, t„, is obtained from the 
decay constant of exponential fits to the same T(t) as the 
tenth row. The twelfth row is the expected universal ratio 
C vv obtained from the previous two rows; note that is is con- 
sistent with being constant except for the largest M data. 
The thirteenth row is the correlation length £' v obtained from 
the velocity autocorrelations £' v = (-j^) 1 ^ , with r v taken from 
the eleventh row. can also be obtained from F(0) and v via 
Eq (|i(i|). This data is presented in row 14. The fifteenth row is 
the force amplitude obtained from Af ~ (F — F^(M)")(£')]/" . 
The sixteenth row is obtained from the finite-size variations 
in the critical force shown in Fig. and the seventeenth row 
is the putative universal ratio Ca inferred from the previous 
two rows; note it is consistent with being constant within the 
(admittedly large error bars). In the eighteenth row, the ve- 
locity amplitude A v is obtained from the fifteenth row and 
from the third row. 



1 


M 


0.0 


0.2 


0.4 


0.6 


2 


A h 


0.98(1) 


1.03(1) 


1.13(1) 


1.23(1) 


3 


A h /A t 


0.39(2) 


0.47(4) 


0.60(7) 


1.0(1) 


4 


At 


2.5(3) 


2.2(4) 


1.9(6) 


1.2(6) 


5 


i 

A> A 13 

K h. 


0.51(6) 


0.6(1) 


0.7(1) 


1.0(2) 


(i 


C K 


0.70(8) 


0.69(8) 


7(1) 


0.6(2) 


7 


A f 


0.16(2) 


0.17(2) 


0.19(4) 


0.21(5) 


8 




1 3(2) 


1.53(3) 


1 8(6) 


3.0(9) 


9 


A v (sim.) 


1.20(1) 


1.49(2) 


2.29(2) 


4.26(3) 


10 


n t 

' V V * X) 


158(5) 


146(10) 


145(16) 


210(20) 


11 




83(1) 


88(2) 


94(4) 


222(6) 


12 


c 


1.90(6) 


1 7(1) 


1 5(2) 


1 1(1) 


13 




10(1) 


11(2) 


13(3) 


31(10) 


14 


tv 


14(2) 


21(3) 


23(4) 


42(11) 


15 


A F 


0.11(2) 


0.10(2) 


0.08(3) 


0.08(5) 


16 


CaA f 


0.21(2) 


0.22(1) 


0.13(1) 


0.14(2) 


17 


C A 


1.9(1) 


2.20(9) 


1.6(2) 


1.8(4) 


18 


A v 


1.7(2) 


2.3(4) 


3(1) 


5(2) 



15 



The correlation time r v can also be obtained from the 
truncated velocity fluctuations, 



Integrating over time yields, 

/ T(t)dt « C vv t v v 2 —^ 



(41) 



(42) 



with C vv a universal coefficient. The resulting t v can, 
together with £„ be used to check other scaling relations. 
See Table I. Note however that there are difficulties as- 
sociated with the subtraction needed to obtain the trun- 
cated correlations in the most interesting regime in which 
the fluctuations are large and an accurate extraction of 
v problematic. 

The velocity-velocity correlations can also be used to 
probe the nature of the dynamics; in particular by study- 
ing the power spectrum which is the Fourier transform, 
r(a>) of T(t). To understand how this is expected to be- 
have in the dissipative limit, it is useful to consider the lo- 
cal velocity-velocity correlation function, (u(x, t)v(0, 0)). 
At long distances, |x| ^> f„ or time separations, t^$>r v , 
this will approach v 2 . But within a correlation space- 
time volume, the local velocities will be characteristic 
of avalanche events and hence be fractal. The correla- 
tions will be proportional to v times a conditional ex- 
pectation of i/(x, t) given that there is motion — i.e. a 
jump — at (0,0). These conditional correlations within 
the space-time correlation volume will reflect the fractal 
structure, being of order l/|x| z ~^ or l/f 1_ ~ whichever 
is smaller. Integrating the associated scaling forms over 
x and Fourier transforming in time, one finds that for 
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with Cps a universal coefficient. 

The log-log plot of the square root of the velocity power 
spectrum (Fig. |l^(b)) appears to exhibit power-law be- 
havior for large lo for both the values of M shown. For 
M = 0, the observed exponent is close to the expected 
value of ~ 0.86 over two and a half decades in fre- 
quency. For M = 0.6, the best fit slope is somewhat 
larger, but some curvature is evident and consistency 
with the dissipative result is not ruled out. In both of 
these sets of data, there is a crossover at low frequen- 
cies to a flat spectrum. This is more pronounced in the 
M = data which are effectively further from the critical 
regime — £/L smaller — than the M = 0.6 data because 
the two sets of data were taken at the approximately the 
same v which is closer to the corresponding v min (L, M) 
for M = 0.6. 

Not only does the value of Fj(M) provide us useful 
information, but its variations do as well. In the finite- 
size limited scaling regime in which L <C £, £ is replaced 
with L in scaling laws and we expect 
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Figure || demonstrates this scaling with the system 
length, the obtained correlation length exponents being 
0.78 ± 0.02, for M = and for M = 0.6, v = 0.81 ± 0.02. 
Calculations of Narayan |l9| yield the leading irrelevant 
eigenvalue at the quasistatic fixed point as approximately 
— |. This suggests a fit of the data of Fig. || with 
the form var(Fj(M)) = CL~i /(l + B v L~i) yielding 
v = 0.77 ± 0.03 with B v = 0.78 ± 0.01. 

We have found that a variety of properties of the steady 
state moving phase as well as the variations of the critical 
force at which the system stops on decreasing the drive 
are all consistent with critical behavior that is indepen- 
dent of the magnitude of the stress overshoots. Given the 
local hysteresis that is intrinsic with stress overshoots, 
this universality is more than a little surprising. In the 
next section, we consider to what extent this applies more 
generally, in particular as far as macroscopic hysteresis. 
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FIG. 10. (a)Time dependent fluctuations of the instan- 
taneous spatially averaged velocity v(t) for M = 0.0, 0.2, 0.4, 
and 0.6. For all data, L = 128 and f™ ax = 0.5. Each sample 
is first equilibrated at some applied force which is then sub- 
sequently lowered at a rate of 4.2 x 10 -8 (see text), until 
the applied force is such that the steady state average veloc- 
ity is approximately 0.04. (b) Log-log plot of the square-root 
of the velocity power spectrum for M = 0.0 (lower plot) and 
M — 0.6 (upper plot). The data have been smoothed by aver- 
aging over groups of five frequencies. The dashed line in both 
plots is the theoretical expectation in the scaling regime: a 
slope of -?±f « -0.86 for M = 0. A fit yields a slope of 
-0.88 ± O.Ollbr M = 0. But for M = 0.6, the slope is rather 
larger than expected: —1.12 ± 0.01. 



V. HYSTERESIS 

We now turn to an analysis of the hysteretic phenom- 
ena that are implied by the coexistence of moving and 
stationary solutions at the same force in the presence 
of stress overshoots. A crucial question which we must 
address is whether hysteresis persists in macroscopic sys- 
tems that are not prepared in special ways. In particular, 
are there hysteresis loops with a width that is non-zero 
in the limit of large systems? If not, as we shall see is 
the case, how does the hysteresis depend on system size? 
Can one understand this in terms of the purely dissipa- 
tive dynamics that appear to control the properties of the 
steady state moving phase? Or is new physics needed? 

In Fig. 0(a), hysteresis loops are shown for typical 
samples of size 128 2 with f™ ax — 0.5 and M from to 
0.6. An upwards arrow indicates the force f J h st , at which 
the system starts moving again after it has been stopped 
at by a gradual decrease in the force. For M = no 
hysteresis is apparent while for positive M, the difference 
~ F c appears to be close to M/4 = M/Z, the 
magnitude of the stress overshoot. On the basis of these 
data, it would appear that the situation is rather simple: 
on decreasing the force the steady state moving phase 
and the stopping process are not qualitatively dependent 
on M, but once the system has stopped, an increase of 
the force by M/Z is required to start it up again. Once 
restarted, the velocity rapidly increases to that of the 
apparently unique moving "state" . The reason for this 
macroscopic hysteresis would appear to be simple: If the 
force on each of the segments caused by the last motion of 
its neighbors before the system stopped was not enough 
to cause it to move, then at later times the force will 
be less than that needed to make a segment move by 
at least M/Z as the stress overshoots from its neighbors 
jumping will no longer be in effect. If this applies to all of 
the segments, it should be necessary to increase the force 
back up again by at least M/Z before anything can start 
to move. This would imply truly macroscopic hysteresis 
that is independent of size for large systems. 



A more careful examination of both the data and the 
argument above shows that it is fallacious: even with 
jrnax _ q a sma ]\ fraction of samples have substan- 
tially narrower hysteresis loops. There must thus be some 
segments that can be restarted by an increase in the force 
from the stopped state by less than M/Z. In the next 
subsection, we discuss the origin of this effect, but first, 
we present and analyze the numerical data. 



A. Distributions of F^ h 8t 



For reasons that will become clear later, we can ob- 
tain more useful data on the hysteresis by increasing the 
strength of the randomness. Most of our detailed hys- 
teresis data is for f™ ax = 1.0 and M = 0.8, the latter 
being sufficiently large that the effects of overshoots are 
strong, but not so large (for this larger value of /™ aa: ) 
that sublattice effects start to play a role. For these pa- 
rameter values, the mean force at which the system stops 
on decreasing F under the procedure discussed in Section 
IV is 



(Fj) « -0.0665 
with the rms variations about this of 



(45) 



^v&r(F}(f™ ax = 1.0, M = 0.8)) w 0.0006 . (46) 
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FIG. 11. (a)Probability distribution of Fj, Prob(Fj hyst ) 
for M — 0.8 and fpm =1.0 for different system sizes, (b) 
Same as (a) but with M — 1.0. The vertical bar on each 
horizontal axis indicates < Fj > +M/Z. 



In Fig. |Tl|(a), the distributions of Fj t t are shown 
for various system sizes. It can be seen that they are 
much broader — by almost two orders of magnitude for 
L = 128 — than the distributions of F*. The shapes of 
the distributions at first appears rather strange: For the 
smaller system sizes, a substantial fraction of the weight 
is in a narrow peak that has similar width to that of 
the distribution of F^, but is shifted up from this by 
an amount fj hy3t — F* ss M/Z. In the largest samples, 
this peak has completely disappeared and we see that 
the width of the hysteresis loop has narrowed consider- 
ably. The narrowing with size of the median width of 
the hysteresis loops is shown in a log-log plot in Figure 
|l2| . As was evident in the shape of the distributions, a 
crossover length of around L = 20 is seen in these data. 
For small sizes the hysteresis loops have width that is 
typically close to M/Z. But for the large sizes, the typi- 
cal width appears to decrease as a power of L: 
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1.15 ±0.03 



1.35 ±0.02. 



(47) 
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for M — 1.0. Note that these exponents are obtained 
from a narrow range of length scales and crossover be- 
havior is likely to be playing a role; we will return to the 
issue of crossovers later. 
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FIG. 12. Log-log plot of the media width of the hysteresis 
loop for M = 0.8 and M = 1.0 with fpm = 1.0. The slope of 
the large L data for M = 0.8 is -0.87±0.02 with an amplitude 
of 3.4 ± 0.2. For the M = 1.0 data, the slope is -0.74 ± 0.01 
with an amplitude of 2.9 ±0.2. We did not include the L = 16 
data points in the linear regression as there is a clear crossover 
at small scales. 



Nevertheless, in spite of uncertainties in the asymp- 
totic size dependence of the hysteresis loops, the overall 
trend is clear: in the limit of large systems, the width of 
the hysteresis loops vanish). This in spite of the fact that 
there are many linearly stable static configurations that 
coexist with the moving state up to the force F c q which 
is substantially greater than FJ h ot . 

Before discussing this result, it is instructive to con- 
sider what happens in the dissipative limit, M = 0. Al- 
though it might appear that there would be no hysteresis 
in this case, this is not strictly correct for finite size sys- 
tems that have been stopped from a moving state. At the 
force at which the system stops, it gets stuck in a some- 
what anomalously strong pinning region (how anomalous 
depends on the rate of decrease of the force) . Getting it 
unstuck from such a region, in the sense that all parts of 
the system move for at least some distance, can require 
an increase in the force that is comparable to the width 
of the distribution of F*. Thus we expect "hysteresis 
loops" in the dissipative limit to have a width of order 

If the systems with stress overshoots behaved like the 
dissipative case in all universal aspects, one would expect 
the asymptotic large system-size dependence of the width 
of the hysteresis loops to have the same exponent as the 
dissipative case, i.e. that [i = v. Up to questions about 
estimation of uncertainties in the presence of complicated 
crossovers, it appears that this is not the case: we seem 
to find 
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corresponding to system size dependence of FJ h st being 

slower than that of F}. If this were indeed the case, we 
would expect that it would most likely hold asymptoti- 
cally for any M > 0. Unfortunately, the range of data 
is not so large as to conclusively rule out equality rather 
than inequality, although if fj, = v one would probably 
need either a large dimensionless amplitude ratio between 
the coefficients of the size dependence of the two critical 
forces, or strongly nonmonotonic behavior; we will later 
explore such scenarios. But for now we focus on the scal- 
ing behavior that seems to be emerging for system sizes 
larger than of order 20 or so for {f™ ax = 1.0, M = 0.8). 

It appears that the larger system sizes do exhibit scal- 
ing behavior of the distributions, a more stringent test 
than exponents. Indeed, the size dependence of the dis- 
tributions provides a useful way to understand the causes 
of the size dependence of the hysteresis. 

Let us assume that restarting on increasing the force 
after stopping occurs via some kind of nucleation process 
whose occurrence is dominated by scales that are much 
smaller than the system size. Then we expect a density 
of nucleation segments — or "seeds" — with a distri- 
bution of values of the local critical forces F s needed to 
restart. If this distribution extends down to F^, larger 
systems are more likely than smaller ones to have a seed 
with a small F s . As the lowest F s in a given stopped 
configuration will be the one that determines F} , this 

lJ Mryst ' 

will yield a distribution of critical forces that becomes 
squeezed down to F* as L — > oo. A simple check on 
the assumption of locality of the seeds is accomplished 
by estimating the distribution of F 'J t for a sample of 
size L, by considering it as being made of b d independent 
samples of size L/b whose Fjs are drawn independently 
from the observed distribution for these smaller size sam- 
ples. In Figs. [l3|(a) and |lj(b), this is carried out for 
L = 128 and b — 2,4, and 8. As can be seen, the distri- 
butions obtained agree very well with those measured for 
the L = 128 samples directly. This agreement is partic- 
ularly striking given that the data for Fj h ^ for L = 16 
in Figs. |ll|(a) and |Tl](b) have a very different form than 
those for the large samples. (This difference in form for 
the distributions is the source of the crossover observed in 
the size dependence of the width of the hysteresis loop.) 
The agreement of the actual distribution of F} for 

G c hyst 

L = 128 with that obtained from the distribution with 
L = 16 suggests that the nucleation process in a sample of 
size 128 is typically dominated by regions whose diameter 
is less than 16. Once such a small nucleation region gets 
going, it will typically expand to make the whole system 
restart, independent of the existence or lack thereof of 
seeds in other regions, or of other stochastic properties 
of the rest of the system. Before analyzing the conse- 
quences of this, we must caution that a crucial question 
is whether such a relation between the distributions of 
Fj b st for systems of size L and L/b holds, in the limit of 



large L, for only a limited range of b, for any b <C L, or 
for b up to some (subdominant) power of L. 

For now, we will extract the shape of the distribution 
from the observation that there is at least a substantial 
range of b over which the relationship between the distri- 
butions of Fj for size L and size L/b does hold. 
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FIG. 13. (a) Distributions, Prob(F] hyst ) , for M = 0.8, 
fpm = 1.0 and L = 128 inferred from data for smaller sizes. 
The L = 128 curve is the same as in Fig. ( [ill ) while the other 
distributions are obtained from the L = 64, 32 and 16 data 
from Fig. ( [Til ) rescaled as described in the text, (b) Same as 
(a) except with M — 1.0. 
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The basic picture of the restarting being controlled by 
the least pinned of many independent seeds enables one 
to relate the size dependence of the median FJ h ^ to the 
form of the distribution. Picking the minimum of the 
FJ h st 's from S = b d subsystems each of which has a dis- 
tribution of Fj that vanishes as (Fj - (F})) dl *- 1 for 

Chyst v Chyst x CI/ 

small Fj h st — (Fj), yields a power law decrease with L 
of the width of the distribution, and, indeed, the actual 
form of the distribution: We expect a Weibull distribu- 
tion with one non-universal scale parameter, fh, E0[. 
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As can be seen in Fig. |l4|, this yields a rather good fit 
to the data for L = 128 with the value of \i extracted 
from the size dependence of the median. If, instead, we 
do a best fit to the shape of the distribution for the 
largest size, we find, /i « 1.08 for M — 0.8 and for 
M = 1.0, f™ ax = 1.0, (j, w 1.25. Note that these values 
are slightly smaller than those obtained from the size de- 
pendence and thus somewhat closer to v; this may well 
be a sign of slow crossover to asymptotic behavior that 
is like the dissipative limit. 
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14. Left plot shows distribution of FX ^ with dark 



solid line a Weibull distribution, Eq. (plf), with scale param- 
eter of 0.048 and fj, = 1.15 obtained from Fig. and dashed 
line a best fit to the Weibull form with two fitting parameters. 
The light solid line is the data for L = 128 and M = 0.8. The 
right plot is, similarly, M = 1.0 with the scale parameter for 
the solid line Weibull distribution 0.078 with [i = 1.35 from 
Fig. n2j, and the dashed line the best fit Weibull distribution. 



B. Origins of seeds for restarting 

We next develop an understanding of the origins of 
the unusual hysteretic behavior that was found in the 
numerical studies: in particular, the origins of the seeds 
for nucleation of motion when the force is increased back 
up again after the motion has stopped. To do this, we 
need to understand how the manifold stops moving as 
F is decreased to F}(M) as it is this that sets up the 
configurations in which the seeds exist. We first analyze 
the basic role of the stress overshoots in the steady state 
moving phase. 

A very crude approximation to the effects of the stress 
overshoots is to ignore their local and transient natures. 
We thus consider an artificial model in which if any seg- 
ment has moved on the previous time step, the force on 
all the segments is increased by M/Z above what it would 
be with purely dissipative dynamics. As long as some- 
thing is always moving, this is identical to merely increas- 
ing the applied force by M/Z. But once the system has 
stopped — because of a decrease in the applied force or 
because of running into a strongly pinned region — no 
segment can move again until the force is increased by 
M/Z; this is because any segment that could move with 
less of an increase in F , should, a fortiori have moved 
already because of the stress overshoot that was present 
before the motion stopped. In terms of the distribution of 
the total force, {/(x)}, on a segment, which must be pos- 
itive for it to jump, the stopped configuration will have 
a depletion layer, no segments with /(x) in the interval 
(-M/Z,0) (and of course none with /(x) > 0). The be- 
havior of this crude model is thus very simple: the v(F) 
curve is shifted down (in F) by M/Z, the steadily moving 
states are identical to those at M = with F — > F+M/Z. 
When the force is increased after stopping, no motion 
will occur until the depletion layer disappears, therefore 
Fj h at = F} + M/Z for each sample. 

Before we return to the model or primary interest, it is 
worth noting that a model which is much less patholog- 
ical than the crude model discussed above nevertheless 
has much of the same behavior. This non-additive stress 
overshoot model has nearest neighbor stress overshoots 
that last for one time step and a "self overshoot" (analo- 
gous to inertia) that likewise lasts for one time step. But 
the stress overshoots are non-additive so that any site has 
a stress overshoot that is either zero or 4r. For example, 
if two nearest neighbors of a segment jump, that segment 
feels a stress overshoot of M- , in contrast to the over- 
shoot it would feel in our primary model. The partial 
equivalence between this non-additive stress overshoot 
model and the crude model can be simply understood: 
If the total force on a segment does not change from the 
previous time step it cannot jump at the next time step. 
Therefore changes in the total force on a segment are 
what determines whether a segment jumps or not at the 
next time step. At fixed F, changes in the total force on a 
segment arise from nearest neighbors jumping at the pre- 
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vious time step and from the segment itself jumping at 
the previous time step. In the non- additive model, each 
of these will involve an extra 4f stress. Thus as long as 
motion has existed somewhere in the systems for more 
than one time step, given a configuration of the crude 
model and which segments have jumped on the previous 
time step, there is an exactly equivalent configuration of 
the non-additive model which will have the same dynam- 
ics at all future times as long as F remains fixed. But the 
dynamics are not fully equivalent: when the first segment 
moves in the crude model, it can trigger others far away; 
this cannot happen in the non-additive model, Neverthe- 
less, the steady state velocity as a function of force will 
be the same in these two models, with the critical force 
in the infinite system limit shifted down by exactly 4f- 
from the dissipative case. The hysteresis loops will also 
be similar, but not identical: in both cases there will be a 
depletion layer of width 4?- after the system has stopped 
and the force will have to be increased by this for motion 
to start again, But in finite size samples the behavior 
will be slightly different as it is much more likely in the 
non-additive model that motion could start in one region 
but die out: the actual critical force for restarting would 
then be slightly higher. The dynamics of the transient 
motion on restarting would also differ due to the locality 
of the non-additive model. 

The decrease of F} with M in the crude and non- 
additive models is the underlying cause for the linear 
decrease of Fg{M) in the primary model: As long as 
segments only move in response to their neighbors mov- 
ing, what is crucial in determining F^ is how the system 
gets through potential sticking points. Some of these are 
likely to involve only one neighbor of a segment moving at 
the previous time step; if they do, then the critical force 
at which they can proceed will, in the absence of other 
changes of the dynamics due to M, be just M/Z lower 
than it would be with M = 0. In the limit of small M, 
we expect that the sequence of jumps at F will be very 
close to that at F + M/Z in the absence of overshoots. 

The non-linear part of the dependence of F} on M 
is of a different origin. As M grows, the equivalence 
between the sequences of jumps at different values of M 
no longer obtains because of, for example, the effects of 
two neighbors jumping at the same time which increases 
the stress on a segment by 2M/Z. This will tend to make 
stopping less likely as a region can be restarted by motion 
in other regions that is caused by such multiple-neighbor 
jumps. As one would thus expect, Fg decreases faster 
than linearly as M is increased. 

The focus on particular sites and whether they can 
be retriggered by a given increase in F is also useful for 
understanding the hysteresis in the model of primary in- 
terest. A crucial question about the local dynamics is: 
How close can a segment be to moving without one of 
its neighbors having moved on the previous time step? 
We must consider the most recent time in the past, say 
time one, at which a neighbor of the segment x of inter- 
est moved. For simplicity, let us assume that none of the 



neighbors moved at time zero. At time zero, the total 
force on x is then 

/(x,t = 0) = F-/ p (0) + <7(0) <0 (52) 

with ct(0) = cr(x,0) and f p (0) = / p (x(0), h(x(0))) the 
initial elastic force and initial pinning strength, respec- 
tively, at x. If n out of the Z neighbors jump at time 
one, 

/(x, 1) = F - / p (0) + a(0) + n(l + M)/Z . (53) 

If this is negative, then x will not move and, at later 
times, the force /(x, t > 1) < —nM/Z as the stress over- 
shoot will no longer apply; thus segment x will not be in 
the depletion layer. If, however, x does jump at time two 
in response to its neighbors jumping, i.e., if /(x, 1) > 0, 
then the total force on it at later times will be 

/(x,2)=F-/ p (l) + £7(0)-l + n/Z (54) 

with a new random pinning force / p (l). As long as 
/(x, 2) < 0, then this segment will not move further 
unless one of its neighbors does. Thus this represents 
a possible local configuration when the system has just 
stopped. 

The condition /(x, 0) < implies, from Eq. ( |52"|) and 
Eq. (H) that /(x, 1) < n(l + M)/Z. The maximum 
of /(x, 2) is then obtained when / p (l) is minimal (i.e. 
zero) and / p (0) maximal (i.e. fp lax ); this yields /(x, 2) < 
fp lax — 1 + n/Z. If the force is now increased back up by 
AF so that /(x, 2) + AF = 0, the segment x will jump 
and could trigger restarting of the overall motion. We 
must thus ask how close to zero /(x, 2) can be. 

If the last motion in any region were always via sin- 
gle neighbors triggering each other, then there would be 
a depletion layer on stopping and macroscopic hysteresis 
for fp 71 ^ < 1 — ^ = | in our case with coordination num- 
ber Z = 4. [This would obtain if the stress overshoots 
were not additive but instead such that any number of 
neighbors jumping yielded the same value of the stress 
overshoot as if only one had.] But in reality, it is possi- 
ble that any number of neighboring segment on a segment 
could jump at one time and then not again; some results 
on the simultaneous hopping of a number of neighbors 
are presented in Table II. As such sets of simultaneous 
jumps can occur for any n < Z — including for n = Z 
even as the system is stopping, there will be no depletion 
layer even for arbitrarily weak randomness and no hys- 
teresis in the infinite system limit. Nevertheless, for weak 
pinning, the depletion layer will only be filled by simulta- 
neous jumps of multiple neighbors followed by a jump of 
the central segment that does not trigger further jumps 
of any of its neighbors, Although we expect that a local 
condition such as this will always occur for some finite 
fraction of the segments, it appears likely that the poten- 
tial seeds for restarting with small increases in F will be 
very rare for weak pinning; this, indeed, turns out to be 
the case. It is the proximate cause of the long crossover 
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lengths apparent in the size-dependent distributions of 

c hyst 

TABLE II. The probability of n = 0, 1,2,3 or 4 of the 
nearest neighbors of a segment having jumped at the previous 
time step given that the segment jumps. For n — 0, there are 
large fluctuations so we record a range of values. The size is 
L = 128 and average velocity approximately 0.1. 
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The above analysis gives a qualitative explanation for 
the lack of macroscopic hysteresis. But to explain the ob- 
served dependence of the widths of the hysteresis loops 
on system size, we must understand the density of states 
of segments with small negative total force on them in 
the stopped configurations, and what happens after mo- 
tion is triggered by one of these seeds as the force is 
increased. Quite generally, the continuous nature of the 
distributions of yield strengths and the discrete nature 
of the stress transfer means that the density of states 
for local properties should either be zero or be strictly 
positive. Specifically, from the above we expect that the 
density of states of the local forces will be positive at 
zero in the stopped state. To check this, we have com- 
puted the probability density per site, r(F), that some 
motion is triggered with a small increase in F to F + dF 
in stopped systems; this is normalized so as to include 
only those samples that have not yet restarted macro- 
scopically (although they might have already had some 
transient local motion). This rate of triggering of jumps 
appears to go to a constant as F decreases to Fj and 
exhibits relatively weak dependence on F as the force is 
increased; we will present the data later in the paper. 

From the data for the distribution of -Fj hyst , which van- 
ishes approximately linearly at F} , a constant density of 
states for triggering motion is perhaps surprising. The 
reason for this must lie not in the seeds themselves, but 
in how they grow. In particular, very close to Fj local 
triggering must be less likely to induce restarting than 
it does at higher forces. In order to understand this, it 
is necessary to investigate the avalanche dynamics, i.e., 
the transient motion in response to triggering of one seg- 
ment. Before considering this in the context of the hys- 
teresis loops and restarting, we analyze avalanches that 
occur in the approach to depinning from below. 



VI. AVALANCHE DYNAMICS 



In the previous section we have seen that to under- 
stand the hysteretic phenomena observed on cycling the 
force up and down, we need to understand how macro- 
scopic motion starts once it has been triggered by a local 
instability that leads to one segment jumping. Before 
studying the case of interest for hysteresis loops, which 
involves initial conditions that are set by the stopping 
process, we analyze the behavior as the force is slowly 
increased from far below the depinning transition start- 
ing from more generic initial conditions. We will call this 
initial depinning. In particular, we are interested in the 
behavior as the depinning transition is approached from 
below. 

Even though there is no steady state motion in this 
regime, there can be local, transient motion in response 
to small increases in F. Such avalanches will not per- 
sist indefinitely for small F because the pinning forces 
in other regions will eventually dominate as long as 
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F < Fj ^M), the — possibly history dependent — crit- 
ical force on increasing F. 

In this section we investigate the dynamics that re- 
sult when F is increased adiabatically: initially by just 
enough that one segment moves. This can then trig- 
ger other segments to hop forward, while F is held fixed 
until the avalanche stops. The same procedure is then 
repeated until F = FJ . For an infinite system, this is 
defined as the force above which the motion persists in- 
definitely in the absence of any further increase. In finite 
systems, there are some ambiguities in how it is defined; 
we choose to define it as the lowest force at which all 
of the segments move during a single avalanche. The 
primary quantities of interest below the depinning are 
the sequence of avalanches and their statistical proper- 
ties: numbers, sizes, durations, etc. More macroscopic 
quantities, such as macroscopic responses to a small but 
non-infinitesimal increase in F can be determined by in- 
tegrating over the properties of the avalanches. 

There are various measures of the size of an avalanche. 
Three of these will be of particular interest. The moment, 
to, of an avalanche is defined as the total motion that 
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this is the quantity of primary interest for earthquakes. 
Alternatively, one can consider the area, a, (in the two 
dimensional case of interest): the total number of seg- 
ments that move at least once during the avalanche. 
Lastly, is the linear size, £, of an avalanche, one mea- 
sure of which is its diameter defined, for example, either 
via some weighted sum of distances of moving segments 
from its center, or as the diameter of the smallest circle 
that will enclose the avalanche. 



A. Scaling 

For the purely dissipative case, M — 0, the scaling of 
the avalanches is related to that of the various quantities 
— h, x, F — Fj and t — discussed in the context of 
the moving phase. In particular, if an avalanche has a 
diameter £, its area will scale as £ df , its duration as £ z , 
the typical maximum displacement — change in h — as 
£'', and its moment as 
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As long as the dimension is less than the upper criti- 
cal dimension, d c = 4 for short range interactions, the 
avalanches will not be fractal and hence 



df 



(57) 



so that the area is a good surrogate, which we will use, 
for the length scale of an avalanche: 



(58) 



Well below the depinning transition, most avalanches 
are small. But as the transition is approached from be- 
low, larger ones become possible, although the distribu- 
tion of their diameters is cutoff by a correlation length £ 
that, in the dissipative limit, diverges at F c q as 
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with the same exponent v as determines the scaling of 
the physically different characteristic length in the mov- 
ing phase, the velocity correlation length. At the crit- 
ical point in the dissipative limit, the distribution of 
avalanche sizes is a power law: 



Prob[diameter > £] ~ — 



(60) 



and a similar relation applies for other measures of size; 
for example, for the area, the exponent is simply changed 
to n/d. Near to the critical force, the distribution of the 
areas of large avalanches has the scaling form 



Probfcfa] = p(a; F)da 



ad. a 



(61) 



where, of the avalanches that occur within a small force 
interval around F, p{a, F)da, is the fraction that of these 
that whose area is between a and a + da [jjj . The scaling 
function V(y — > oo) decays rapidly while for y — > 0, it 
goes to a constant. 

The statistical "tilt" symmetry of the system that was 
used earlier to yield the scaling law, i = 2 — £, can also 
be used, via relating the polarizability to the avalanche 
production rate and the distribution of their sizes, to 
show that for M — 0, 



n=d--=d-2+( 
v 



(62) 



as derived in Ref. 0. The one crucial assumption is 
that the rate, r(F), of avalanche production, defined as 
l/AF times the number of avalanches per unit area of 
the system as the force is increased by a small amount 
from F to F + AF, tends to a finite non-zero constant 
at the critical force. 

We now turn to an analysis of the data for avalanche 
statistics and properties of the avalanches: for the dissi- 
pationless limit, to check the theoretical predictions out- 
lined above; and for non-zero M, to investigate the effects 
of stress overshoots. 

An easy quantity to measure is the cumulative distribu- 
tion of all the avalanches as the applied force is increased 
to Fj AM). This is given by 
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Assuming that r(F) approaches a constant as F — > , 
the scaling laws for the dissipationless case yield a power 
law cumulative avalanche area distribution with an ex- 
ponent of unity, 



K cum (M = 0) 



1 



(64) 



independent of the values of the other exponents. This 
thus provides a good test of the general scaling theory 
that does not depend on particular predictions for expo- 
nents. Although we do not expect the universal aspects 
of the avalanche statistics to depend on details of the 
initial conditions, to avoid effects that might arise from 
smoothening out rough initial conditions, we take the ini- 
tial configuration to be approximately flat: specifically, 
the initial {/i(x)} uniformly distributed in the interval 
[0,1]. 



B. Dissipative limit 

We first analyze the data for M — 0.0. In Fig. |l5|(a), 
the cumulative avalanche area statistics are shown; a fit 
to the data on a log-log plot yields an exponent some- 
what less than the theoretical expectation of K cum = 1 , 
If we restrict consideration to those avalanches that occur 
in the region close to the critical force in which most of 
the activity occurs, specifically, within the applied force 
region of [FJ — 0.1, ], the apparent exponent is 
roughly the same, K cum = 0.85 ± 0.03 as shown in Fig. 
|l5|(a) . Before trying to understand the apparent discrep- 
ancy of this with the scaling prediction, we consider the 
statistics of avalanches that occur in the critical regime, 
specifically, only those that occur for F > (F+(M)). The 
distribution of these also decays as a power of the area, 
as shown in Fig. p"5|(b). But the power is much smaller: 
0.37 ± 0.02. Because the correlation length that would 
cutoff the avalanche distribution is of order the system 
size in this regime, the distribution should essentially be 
that of critical avalanches with an exponent n/d. We 
thus obtain an estimate 



k(M = 0)w 0.74 ±0.04 



(65) 



to be compared with the theoretical expectation of n = 
d — 2 + £ = d — i equal to £ in two dimensions. We 
see that the agreement with our data for £ is quite good, 
suggesting that the basic scaling scenario is correct. 

There are various possible sources for the substantial 
discrepancy of the apparent cumulative exponent from 
unity. As can be seen in Fig. |l5|(a),(b), the rate, r(F) of 
avalanche production for M = increases sharply as the 
critical force is approached. Although it does appear to 
go to a finite constant, as it should on general grounds 
since r(F) is a locally determined property, the precursor 
increase will bias the cumulative statistics as there are 
more avalanches produced near the critical point than 
further away, and these are the ones that have possibili- 
ties of being large. This will tend to put more weight in 



the large avalanche part of the cumulative distribution 
thereby decreasing the apparent K cum . [The extreme 
limit of this weighting of those near the critical force 
would just yield the critical exponent n/d instead as in 
Fig. n3(b).] In particular, if there is a cusp singularity in 
the rate of avalanche production as F — > FJ with an 
exponent a, then this would induce a multiplicative cor- 
rection to the cumulative avalanche size distribution of 
the form l — C/a~ which could complicate interpretation 
of data. 

In general we expect there to be corrections to scal- 
ing arising from weakly irrelevant operators at the RG 
fixed point that governs the depinning critical behav- 
ior. In terms of areas, the leading irrelevant eigen- 
value —0m |l{| would give corrections of the form 
a~ Kcum (l + BaT e l d ) with 6/d w | in two dimensions. 
Such corrections would dominate over those from a cusp 
in the avalanche production rate unless a < 0.5. Indeed 
the simplest expectation is that the cusp in the avalanche 
production rate is controlled by exactly this correction 
exponent, 
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FIG. 15. (a)Log-log plot of the probability of avalanche 
area a occurring within the interval [a/y2, y2a) as F is in- 
creased to Fl awl {M) for f™ ax = 0.5. The circles are for 
M = 0.0, squares M = 0.2, triangles M = 0.4, and inverted 
triangles M = 0.6. The open symbols represent L — 128 
data, while the closed symbols represent L = 64. The asymp- 
totic slope should be —K cum /d; fits yield: for M = 0.0, 
-0.89 ± 0.02; for M = 0.2 , -1.13 ± 0.02; for M = 0.4, 
-1.25 ± 0.01; and for M = 0.6, -1.53 ± 0.02. The dashed 
curve is the result of a two-parameter fit to the M — data: 
Ca- 1 (l + B cum a- 1/3 ) with B cum = -1.21 ±0.11; with K cum 
also a fitting parameter: Ca~ Kaum (l+Bcumfi' 1 ' 3 ), these data 
yield K cum = 0.83±0.02 and B cum = 0.53±0.01. (b) Same as 
in (a) and but including only the avalanches that were initi- 
ated with an applied force greater than i 7 ^ (M) — well within 
the critical regime. The expected slope is —n/d; for M = 0.0, 
a best fit yields —0.37 ±0.02. For non-zero M, large crossover 
effects are evident. 



"(F) ~ r(F c0 ) ~ C r (F c0 - F) 6 



(66) 



corresponding to a w 0.5. A fit to this form with r(F c o) 
and F c q fixed, and C' r and 9v as the free parameters, 
yields a lower value of a. If we take into account a correc- 
tion to scaling by fitting the log- log plot of the cumulative 
avalanche area distribution with a form that includes the 
leading correction (shown in Fig. |l5|(a)), the cumulative 
avalanche size exponent is changed from 0.89 ± 0.02 to 
0.83 ± 0.02, i.e. in the wrong direction. But if the data 
are fit with the expected K cum — 1 and an a -1 / 3 cor- 
rection, — a fit with the same number of parameters as 
an undetermined power law with no correction — the in- 
ferred correction to scaling amplitude is B w —1.2 ±0.1, 
not unreasonably large. If this were the actual form of 
the distribution, some downward curvature at the largest 
sizes would be expected as seen in the dashed fitting line 
in Fig. 15(a). A competing tendency, however, is the 



system area; this is evident in Fig. [l5|(a) in which data 
for L = 64 and L — 128 are shown. These effects combine 
to make the real uncertainties in the exponents here, and 
probably for other quantities, substantially larger than 
the apparent uncertainties. 

Note, however, that the useful range of length scales 
available for avalanche data and other quantities for 
which one is trying to extract infinite system results from 
finite system data are, because of the crossover when the 
length scales approach the system size, less than is avail- 
able for quantities, such as the variance of the critical 
forces, that are intrinsically properties of finite size sys- 
tems. Thus we might hope to have more confidence in 
exponent estimates extracted from such intrinsic finite 
size properties. 
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FIG. 16. (a)For initial depinning, log-log plot of the 
probability of an avalanche area a occurring in the interval 
[a/V2,V2a) as F is increased to F^^JM) for f™ ax = 1.0 
and L — 128. The open circles are for M = 0.0, the open 
triangles denote M = 0.4, the stars represent M = 0.8, and 
the M = 1.0 data are denoted by x's. The fitted slope for 
M = 0.0 is -0.90 ± 0.01; for M = 0.4, -1.02 ± 0.02; for 
M = 0.8 , -1.08 ± 0.02; and for M = 1.0, -1.14 ± 0.02. 
The dashed curve is a fit to the M = data that includes 
corrections to scaling as in Fig. |l5|(a); yields a correction 
to scaling coefficient of B CU m = —0.97 ± 0.01. (b) Same as 
(a) but including only the avalanches that were initiated with 
an applied force greater than F}(M), in the critical region. 
Comparisons with data taken on stopped samples (starting 
with F > Fj and lowering F until the manifold stops before 
increasing again) . For the initial depinning, the same symbols 
as in (a) are used. For the stopped samples, the right trian- 
gles are for M = 0.8 and the left triangles for M = 1.0. For 
M — 0, the slope is expected to be —n/d — — Q/d w — i a fit 
yields -0.39 ± 0.02. For M > 0, the data do not fit this and 
crossover effects are evident. For comparison, dashed lines 
are shown with a slope of —0.9, which corresponds to the 
measured cumulative avalanche area exponent K cum for the 
M = 0.0 case. 



better and yields a correction to scaling amplitude of or- 
der unity, more precisely —0.91. Note, however, that, 
as for the distribution of avalanche sizes, such fits suffer 
from finite size effects. (3) 
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1. Durations, moments, and roughness 

From the avalanche data, specifically for the durations 
of avalanches, one can determine the dynamic exponent 
z: with non-fractal avalanches, the duration r will scale 
as a'S . In Fig. [l7](a) we plot the mean duration of 
avalanches whose area is within a factor of y/2 of a as 
function of a. The slope of the log-log plot yields 

z(M = 0) ps 1.40 ± 0.02 . (67) 

As for the avalanche statistics, we can attempt to take 
into account corrections to scaling by fitting to the form 

(r(a)) ps i + g a ^- 1/3 ■ This three parameter fit yields z — 
1.52 ± 0.02 and B z = -0.58 ± 0.05; this value of z is 
closer to the value of from a naive extrapolation of 
the 4 — e expansion, although z is likely to have higher 
order corrections in e even if Q does not. Ref. p5[ shows 
that z = 2- |e-0.04321e 2 ps 1.38. The z extracted from 
the avalanche dynamics agrees quite well with this value 
without any corrections to scaling. 

The roughness exponent Q can be extracted from the 
moment as a function of the area: m ~ a~J~ . Factor-of- 
two logarithmically binned histograms are shown in Fig. 
^T| ; for the dissipative case these appear to yield a rather 
small value of 

( a val(M = 0) ps 0.44 ±0.02 (68) 

however there is definitely upward curvature observable 
at the larger sizes and a fit with C = | and an a~ 3 correc- 
tion to scaling (same number of parameters) is somewhat 
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FIG. 17. (a) Log- log plot for the mean duration of 
avalanches whose area is in the interval [a/y/2, y/2a). All data 
are for f™ ax = 0.5, L = 128. The M > curves have been 
shifted along vertical axis for clarity by multiplying by suc- 
cessive factors of 3/2. The asymptotic slopes should be z/d. 
Linear fits for both M = 0.0 and M = 0.2 (not shown) yield 
z = 1.40±0.02 .The dashed line is a non-linear fit that includes 
corrections to scaling; for M = 0.0 it yields z — 1.52 ± 0.02 
with correction amplitude B z = —0.58 ± 0.05;; see text. Lin- 
ear fits for the M = 0.4, M = 0.6, M = 0.8, M = 1.0 data 
yield values of z/2 of 0.69 ± 0.02, 0.70 ± 0.03, 0.69 ± 0.01 
and 0.68 ± 0.02 respectively, (b) Log-log plot for the median 
duration avalanches with area in the interval [a/y/2, y/2a). 



A roughness exponent can be more directly extracted 
from the statistics of the roughness of the manifold cx- 
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actly at the critical force, , In Figs. |l9|(a)(b), fits to 
the spatial power spectrum, (\h(k)\ 2 ), yield C,(M = 0) « 
0.65±0.02 and 0.69±0.02 for /™ aa; = 0.5 and /™ a:E = 1.0 
respectively. These values are close to those observed for 
systems that have stopped after a decrease of F from the 
moving phase. 

The apparent discrepancy between the scaling of the 
typical displacement of avalanches and the roughness at 
the critical point is somewhat troubling, although, as we 
have seen, it can readily be accounted for by corrections 
to scaling. This sort of discrepancy has been seen previ- 
ously in simulations of manifold depinning pl| . In that 
case, the roughness exponent at the critical force was 
similar to ours, but the displacement of the avalanches 
scaled with a larger exponent. The fact that the discrep- 
ancy can be in either direction supports the belief that 
it is due to corrections to scaling rather than a difference 
between the two exponents, as had been conjectured in 



M=0.0, 0.2, 0.4, 0.6 




a 

FIG. 18. Log-log plot of the mean moment of avalanches 
with area in the interval [a/v2, ■\/2a) for L = 128 and 
= 0.5 using same symbols as in Fig. |l|(a). The ex- 
pected slope is (d + Q/d. A linear fit (solid line) for M — 0.0 
yields 1.22 ± —0.02 For other M, the same exponent fits to 
within one standard deviation. The dashed line is a fit in- 
eluding corrections to scaling: Ca^~ (1 + A^a 3 ) with C 
and Bq = —0.91 ± 0.02 fitting parameters. 



v{M = 0) w 0.82 ± 0.04 (69) 

for f r p nax = 0.5 and v{M = 0) = 0.78 ± 0.02 for the 
stronger pinning. The scaling law £ + £ = 2 is consistent 
with the inferred exponents at the level of a few times 
the apparent error bars. 

In spite of the discrepancies noted above for the typical 
displacement of avalanches and their cumulative statis- 
tics, overall the data for avalanches in the dissipative 
limit appear to be consistent with those obtained from 
the moving phase and with theoretical predictions. In 
particular, if one includes corrections to scaling of the 
anticipated form, all the data is consistent. 

C. Avalanches with stress overshoots 

We now turn to avalanche properties for M > 0. Be- 
cause the distribution of avalanche sizes will turn out to 
raise the most questions, we first study the duration ver- 
sus area and moment versus area. 



(a) 




0.1 1.0 10.0 

l*l 



As was done for the critical force defined from the mov- 
ing phase, we can also extract the correlation length ex- 
ponent from the finite-size scaling of the variance of the 
critical force. In Fig. 20, the data are shown and expo- 
nent estimates of v ^ =a) » 1.22 ± 0.04 and 1.28 ± 0.02 
extracted from the data for, respectively, fp ax = 0.5 and 
f™ ax = 1.0. These yield 
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FIG. 19. (a)Log-log plot ol the square root ol the power 
spectrum contrasting the system just below the critical force 
on initial depinning, Fj ,,(open symbols); with its behavior 
when it has "just started" ummediately after every segment 
has moved during a single avalanche event at F} (solid sym- 
bols). Each set of data represents an average over 1000 sam- 
ples, and the curves have been shifted along the vertical axis 
for clarity. The stars and x's denote f™ ax = 1.0, the other 
data f™ ax = 0.5. For the data just below Fj aval , the fitted 
roughness exponents, £, are indicated by the uppermost solid 
lines: with f™ ax = 0.5: 0.71 ± 0.05 for M = 0.0, 0.55 ± 0.04 
for M = 0.2, and 0.02 ± 0.02 for M = 0.6; with f™ ax = 1: 
0.57 ± 0.04 for M — 0.8 . The lower solid curves are fits with 
corrections to scaling of the form cik~^ 1+ ^ (l + B(k 2 / 3 ), with 
fitting parameters ci and B c = -0.006 ± 0.002 for M = 0.0, 
0.034 ± 0.007 for M = 0.2, and 0.283 ± 0.009 for M = 0.8. 
For the just started data, with f™ ax = 0.5, inferred £'s are 
0.65±0.02 for M = 0.0, 0.65±0.02 for M = 0.2, and 0.69±0.02 
for M = 0.5 (not shown); and 0.72 ± 0.03 for f™ ax = 1.0, 
M = 0.8. No fit was attempted for M = 0.6, f™ ax = 0.5. 
The error bars are achieved in the same way as described in 
Fig. |^(a) and the data has been thinned out as in Fig. ^ (a) 
as well, (b) Same as (a) for comparison of amplitudes: the 
curves have not been shifted and a narrower range of |fc| is 
shown. 



1. Durations and moments 

In Figs. 0(a) and (b), the mean and median durations 
of avalanches whose area is within a factor of y2 of a are 
plotted for various M, Except possibly at the largest 
sizes, the data look remarkably similar to the data for 
M = and an exponent of z w 1.4 would be inferred 
from each. In Fig. ^l] the mean moments are plotted; 
again there seems to be remarkably little dependence on 
M . For both the duration and the moment data, not 



only do the exponents look similar to those at M = 0, 
but the amplitudes do also. Given the relatively strong 
dependence of the velocity amplitude on M (Sec. IV), 
this is particularly surprising. 
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FIG. 20. Log-log plot of the rms variations of Fj aval (M) 
as a function of the system length L. Filled symbols denote 
f™ ax = 0.5 and open symbols f™ ax = 1.0. The stars repre- 
sent M = 0.8, f™ ax = 1.0. Other symbols are as in previous 
figures. The slopes of the two M = 0.0 curves are — 1.22±0.04 
(filled) and —1.28 ± 0.02 (open), shown by solid lines. For 
M > 0, crossover is evident and he dotted lines are merely 
guides for eye. 



An amplitude that would be expected to be propor- 
tional to the time- length scaling factor, At, can be ex- 
tracted from the mean duration of avalanches as a func- 
tion of their area. Although these would be expected 
to differ by an unknown numerical factor from those ex- 
tracted from measurements in the moving phase, the ra- 
tios might have been expected to be universal. The am- 
plitude Af al « 1.33 for M < 0.4 and Af val « 1.44 for 



1 aval 
H 

M = 0.6, f™ ax = 0.5. It is apparent that these ampli- 
tudes exhibit substantially less dependence on parame- 
ters than do the A t inferred from the moving phase. Al- 
though at first this seems troubling, it should, perhaps, 
not be: Even in the absence of overshoots, the properties 
of avalanches can depend on how the system is prepared. 
Nevertheless, the differences here are probably related to 
some of the more puzzling differences that occur in dis- 
tributions of avalanche areas and distributions of critical 
forces; these we discuss further below. 
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FIG. 21. Log-log plot of the time evolution of the av- 
erage moment of avalanches initiated between F^(M) and 
Fj AM) (within a small force interval). The normaliza- 
tion factor is the number of avalanches that still moving at 
the time t. The light solid line displaced slightly below the 
M — 0.0 curve is the result of a least squares fit to a power 
law whose exponent yielding an exponent of 1.63±0.01 in the 
scaling regime. In the quasistatic scaling regime, the expo- 
nent should be (£ + d)/z. For all data, L = 128. For clarity, 
< Moment(f) > for M = 0.2 has been multiplied by 2; for 
M — 0.6, the multiplicative factor is 4, and for the M — 0.8 
data, it is 6. 



2. Distributions of avalanche area 

In contrast to the properties of avalanches of a given 
area, the distributions of avalanche sizes depend strongly 
on M. In Figs. [l5] (a) and [l6| (a) the cumulative distri- 
butions of all avalanches are shown for f™ ax = 0.5 and 
ymaz — i q respectively: the data are straight on a log- 
log plot over a substantial range of areas but the slope 
appears to vary continuously with M, for f™ ax = 0.5 
from 0.89±0.02 for M = to 1.53±0.02for M = 0.6 and 
for f™ ax = 1.0, from 0.90±0.01 for M = to 1.14 ±0.02 
for M = 1.0. At large areas, of order 10% or so of the 
system area, the distributions of avalanches falls off in- 
creasingly more rapidly as M increases. Whether this is 
an intrinsic effect or a finite size effect is an important 
issue to understand. 

Data taken for stronger pinning, f™ ax — 1.0, are 
shown in Fig. [l6](a). Surprisingly, these cumulative data 
depend far less on M over the range shown. Although 
the apparent exponent varies slightly, all are probably 
consistent with K cum = 1 with corrections to scaling. 
More strikingly, there is no dropoff seen for large sizes 
even with M = 1.0. 



The dependence of the avalanche distributions on 
stress overshoot are more marked if we restrict consid- 
eration to avalanches that occur for forces above F^, 
and 
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As mentioned earlier, for M — these 
yield an exponent, expected to be equal to n/d = £/2, 
of 0.38 ± .03. But as M increases the distributions ap- 
pear to follow roughly the dissipative behavior for small 
avalanches only to fall substantially below it for large 
ones. This difference is most pronounced for large M: 
the largest avalanches are about an order of magnitude 
rarer for M — 0.4 than for M — for the weaker pin- 
ning and similarly for M = 1.0 versus M = for the 
stronger pinning. In contrast to the data for the cumu- 
lative numbers of avalanches, these data are suggestive 
of typical crossover behavior: from dissipative behavior 
for small sizes to something else for large sizes, with a 
crossover length that is long for small M and shrinks as 
M grows. Note that this form of M dependence of a 
crossover length is what would be expected if M were 
a relevant perturbation about the dissipative depinning 
fixed point. But we must exercise caution before drawing 
such a conclusion. 

From the data of Fig. [TJ| for f™ ax = 0.5, M = 0.6 
it appears that there may be some interesting size de- 
pendence to the crossover: these data for L = 64 fall 
somewhat below those for L — 128 for large avalanches. 
One can test whether the crossover is a finite size effect 
(in contrast to one arising from putative relevance of M) 
by plotting the data versus the scaled area, a/L 2 : this 
is done in Fig. |22|(a) for a wide range of system sizes 
with M = 0.6. It can be seen that the rollover at large 
areas does not scale simply with system size. But neither 
is it consistent with a system-size independent crossover: 
there is much less than the expected factor of 64 differ- 
ence in the crossover value of a/L 2 between L = 16 and 
L = 64. A plausible intermediate conjecture would be 
that the crossover scales as a power 1 — T of system area: 
this could arise from dangerous irrelevancy of M. The 
plot in Fig. ||^(b) versus a/L appears roughly consistent 
with such a conjecture with T w 0.5. But this requires 
the introduction of an extra scale and exponent into the 
interpretation of the data. Another interpretation, one 
that does not require such additional hypotheses, is that 
systems with f™ ax = 0.5, M = 0, 6 are in the midst 
of some type of crossover for the range of system sizes 
investigated. As discussed later, we believe that this is 
probably the case. 



3. Roughness from avalanches 

The absence of large avalanches in the presence of sub- 
stantial stress overshoots for f™ ax = 0.5 suggests that, 
right up to FJ , the manifold should be less rough on 
long length scales than in the dissipative limit. This is 
indeed found to be the case: In Figures |l^(a) and (b), 
the spatial power spectrum of the roughness is shown as 
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a function of wavevector just below Fj . In the dissi- 
pative limit, a power law that is consistent with a rough- 
ness exponent of ( — 0.65 ± 0.02 for f™ ax = 0, 5 and 
C = 0.71 ±0.05 for f™ ax = 1.0 obtains down to the small- 
est wavevectors. In contrast, for f™ ax = 0.5, the spectra 
with stress overshoots are found to become flat at long 
length scales; this is particularly pronounced at the larger 
M. Although we have not studied the system size depen- 
dence of these spectra in detail, for M = 0.6 the rough 
magnitude of the wavevector at which the crossover from 
power-law rough to flat occurs — as seen in the blowups 
of the small |k| regime in Fig. |l9|(b) — is similar to the 
wavevector inferred from the inverse square-root of the 
system area at which the drop-off in the avalanche area 
distributions occurs. But to again confuse the interpreta- 
tion, different behavior is once again seen for the stronger 
pinning, fp Hax = 1, M = 0.8, samples: these show ap- 
parent roughness exponents that are only slightly smaller 
than in the absence of overshoots and the roughness has 
less of a tendency to saturate at small wavevectors than 
the weaker pinning data. 

Information can also be garnered from the system size 
dependence of the average maximum width, IFmax =< 
max(|/i(x) — h(x))\ >, of the manifold at Fj . In the dis- 
sipative case, W max ~ as expected for simple scaling. 
The stronger pinning samples with f™ ax = 1.0, M = 0.8 
show behavior for H^max that is relatively similar to the 
dissipative case, consistent with the spatial power spectra 
(\h(k) 2 \). For the samples with f™ ax = 0.5, M = 0.6, 
however, much weaker size dependence is observed: Fig. 
Bp . We consider several possible reasons that this pecu- 
liar behavior for these parameter values might occur. 

If, as appears to be the case for samples with f™ ax = 
0.5, M = 0.6, the avalanche production were sharply 
cutoff at a diameter of order L 1_T from its distribution 
in the dissipative case, the largest avalanches that occur 
would change h locally by of order L^ 1_T ^, yielding a 
Wmax of this order. If, instead, the crossover were to 
a distribution decaying as a larger power-law of the area 
for large avalanches, the apparent roughness exponent ex- 
tracted from PFmax would be somewhat larger. But if the 
largest avalanches were qualitatively different from dis- 
sipative ones — - perhaps like failed nucleation bubbles 
with W ~ L — this could lead to an increased appar- 
ent roughness exponent for W max . Unfortunately, at this 
point, it is hard to conclude much from this set of data 
except that something peculiar seems to be going on for 
jmax _ q 5^ M — q g g ut peculiarity in this range of 
parameters appears also in other quantities. 
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FIG. 22. (a) Finite-size scaling plot for cumulative 
avalanche area distribution during initial depinning. (b) Scal- 
ing plot of the same data with T = 1/2. 
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5. Variations of FJ 




FIG. 23. Log-log plot of the maximum width of the man- 
ifold just below the initial depinning at Fj . Data is aver- 
aged over 10 samples. For the M = 0.0 curve, the slope is 
0.70 ±0.03; for the M = 0.2 curve, it is 0.56 = ±0.02. For the 
M = 0.6 curve, the result of the least squares fit is a slope of 
0.29 ± 0.06. And finally, for the M = 0.8 curve, the slope is 
0.64 ±0.06. 



4- History dependence of roughness 

The data at Fj for both the spatial power spectrum 
of the roughness and for W max should be contrasted with 
those in the steady state moving phase , Fig. 0(a), and 
with those in the "just started" moving state that has not 
reached steady state: Fig. |n](a). The latter two both 
show power-law roughness with exponents that are con- 
sistent with being universal independent of M and inde- 
pendent of whether the manifold has been slowly stopped 
from the moving phase or has barely started moving on 
initial depinning — albeit with a velocity dependent cut- 
off at long scales in the latter case. In contrast, the data 
just below FJ exhibit roughness that depends substan- 
tially on M. 

The trends with changing M are also different for the 
stopped samples (at Fj) and the samples at initial de- 
pinning (FJ ): When the motion has been stopped 
from the moving phase, the roughness is slightly larger 
for larger M at least for the weaker pinning. In con- 
trast, the long length scale roughness in systems that 
started flat and had F increased to Fj tends to de- 
crease with increasing M because of the rarity of large 
system-roughening avalanches with this history. 



We have seen that there appear to be interesting dif- 
ferences between the physics of avalanches in the dis- 
sipative limit and those with stress overshoots. In the 
dissipative case, the critical force, FJ , defined from 
the first system-sized avalanche, and that, FJ, defined 
by when steady-state motion ceases, are essentially the 
same; the latter being slightly larger and with a some- 
what smaller variance because of the way it is defined 
(the system stops in a somewhat anomalously strongly 
pinned region). See Table III. With substantial stress 
overshoots, in contrast, while the lowest force at which 
a finite size system can have a system-sized avalanche is 
roughly Fj(M), the distribution of Fj extends far fur- 
ther up, into the region in which, once the motion starts, 
it is extremely unlikely to stop again. From Fig. |2^. it 
is seen that for f™ ax = 0.5 and L = 128, the rms varia- 
tions of FJ are about a factor of four larger for M = 0.4 
than those for M — 0, and for fp lax — 1.0, they are a 
factor of ten larger for M = 0.8 than for M = 0, even 
though the corresponding rms variations of Fj hardly dif- 
fer. This factor of ten difference in F = Fj corresponds 
crudely to a factor of 10" w 6 difference in length scale 
suggesting that some scale of order a sixth of the system 
size might appear for the initial depinning history with 
jmax = i o, M = 0.8. Recall that such a scale was also 
apparent in the distribution of the force, FJ h st , at which 
the system restarts after being stopped from the moving 
phase. 

TABLE III. Comparison of Fj(M) and Fj ava[ {M) for 



L = 128. 


M f™ ax 


(Fj) 


(Fc n „„., } 


(var^M))) 1 / 2 


(var(Fj_,(M))) 


0.0 0.5 


-0.0779 


-0.0785 


0.0004 


0.0006 


0.1 0.5 


-0.1087 


-0.106 


0.0004 


0.003 


0.2 0.5 


-0.1404 


-0.137 


0.0004 


0.003 


0.4 0.5 


-0.2078 


-0.199 


0.0005 


0.002 


0.6 0.5 


-0.2866 


-0.267 


0.0004 


0.005 


0.8 1.0 


-0.0665 


-0.060 


0.0006 


0.009 


1.0 1.0 


-0.1513 


-0.141 


0.0006 


0.007 



1/2 



31 



The data for the rms variations of FJ are shown 



in Fig. |20| and histograms of the distributions in Figs. 
U(a),(b),(c). Note the strong deviations from simple 
power law behavior of the intermediate M data for 
the variances. Surprisingly, the large M samples with 
jmax = 10 have size dependence of the variance of Fj^^ 
that is, except for the overall amplitude, quite like that 
of the dissipative limit. This again recalls the puzzling 
behavior that was seen for the distribution of avalanche 
areas: the strong pinning large M systems appear more 
similar to the dissipative limit than do intermediate M 
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FIG. 24. Distributions of FJ , with, from left to right in 
each plot, L = 128, 64, 32, 16. (a) M = 0.2, f™ ax = 0.5. (b) 
M = 0.6, f; nax = 0.5.(c)M = 0.8, f™ ax = 1.0. 



6. Crossover behavior 

The data we have presented on the statistics of 
avalanches and the roughness they induce show compli- 
cated dependence on the system size, on the magnitude 
of the overshoots, and on the strength of the random 
pinning. But collectively they do suggest a plausible sce- 
nario involving a subtle crossover as a function of length 
scale: For small M one would expect dissipative-like be- 
havior for small samples as the effects of small overshoots 
will need a substantial range of length scales to build up. 
At the other extreme, the data for large M and strong 
pinning suggest that very large samples show behavior 
characterized by the same exponents as the dissipative 
limit, but with very different amplitudes. In between, at 
an M dependent scale, there must then be a crossover 
from one to the other amplitude. On, for example, a 
log-log plot of the rms variations of Fj versus system 
size this would show up as a regime of slope ^ , followed 
by a crossover regime of lower slope to an asymptotic 
regime of slope again -. The intermediate M samples 
could then appear to have substantially different scaling 
by virtue of being in the crossover regime over the range 
of length scales studied, In the next subsection we see 
how such crossover behavior might arise from the pro- 
duction and development of avalanches and the effects 
each has on subsequent ones. 
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7. Avalanche production and development 

We have seen that with stress overshoots, static con- 
figurations exist in finite size samples at forces that are 
substantially above (on the scale of variations of F*(L)) 
the force, F*, at which the whole system can move. This 
suggests that in this regime there may be a change in 
character of the avalanches as they become large. In 
contrast to below F} where all avalanches will stop, if 
an avalanche above F* becomes larger than some char- 
acteristic size, naively, perhaps the velocity correlation 
length, Cv(F), in the moving phase, it will runaway and 
the whole system will move. This will strongly affect 
the distribution of avalanche sizes in this regime, yet it 
might not have much effect on the durations or moments 
of those that do not runaway: the evidence discussed 
above is that the durations and moments of these are 
indeed rather similar to the large ones that occur with 
M = 0. 

The rates, r(F), of avalanche production as F is in- 
creased adiabatically are shown in Figs. |24|(a),(b) and (c) 
for a variety of values of M and f™ ax ' . These are normal- 
ized by the number of samples that are still pinned, i.e. 
still below their F 'J at the given F. In contrast to the 
dissipative limit for which r(F) increases sharply with F 
as F — > Fj w Fg , the data for non-zero M show a peak 
in the avalanche production rate somewhat below F} and 
then a marked decrease as F increases into the region in 
which some samples have self sustaining avalanches and 
become depinned. For the strongly pinned systems this 
suppression is strong; the behavior in the region in which 
r{F) is small is shown blownup in the figure. Note that 
in this regime the avalanche production rate is two orders 
of magnitude lower for M = 0.8 and M — 1.0 than for 
M = 0. 

Although we have not studied it in detail, this decrease 
in avalanche production rate is presumably related to the 
depletion of close-to-unstable sites in regions in which 
moderately sized avalanches have already occurred; this 
is analogous to the depletion that occurs after the whole 
system is stopped from the moving phase. Such depletion 
will suppress both the density of "seeds" of avalanches 
and the probability of them becoming large. Most new 
avalanches would be expected to occur in regions that 
have not yet had substantial sized avalanches. The ex- 
ception to this are avalanches that get up enough steam 
that they can run through the depleted regions. Once 
an avalanche does this, it will continue growing until it 
sweeps through the whole system and becomes self sus- 
taining. As the chances of this happening will depend on 
the properties of various regions of the system, how large 
it has to become to runaway is likely to be subtle. 
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FIG. 25. (a) Avalanche production rate during initial de- 
pinning. From top to bottom of the left plot, f™ ax = 0.5, are 
M = 0.0, M = 0.2, M = 0.4 and M = 0.6. The solid line fit 
to the M — data is 2.35 - 3.3(-0.0785 - F) 0/M . The dotted 
line is a fit in which only CV (from Eq. (66)) is allowed to vary 
but a is fixed at 0.5; note there is only one fitting parameter 
for the latter. In the right plot, f™ ax = 1.0, top to bottom: 
M = 0.0, M = 0.4, M = 0.8, M = 1.0. (b) Blow-up of (a) 
near the critical region for f™ ax = 0.5 data. Up arrows along 
the horizontal axis indicate the median Fj^^ , while down ar- 
rows indicate (F}). (c) Blow-up of (a) for f™ ax =1.0 data, 
(d) Comparison of r(F) for initial depinning (insets) with that 
in the hysteresis loop. Note that for samples with L — 128 
(shown) the hysteresis loops are considerably narrower than 
their naive width M/Z. 
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There are three possible behaviors for the runaway pro- 
cess of large avalanches in the regime above Fj : (i) that 
there is a finite sample-size-independent length scale, 
£n(M) above which avalanches typically run away, (ii) 
that the scale at which runaway occurs grows with sys- 
tem size, but as a decreasing fraction of system size, or 
(iii) that the runaway scale is proportional to system size 
in large samples. In the dissipative limit, the last of these 
obtains with a proportionality constant close to unity. 
In contrast, in the artificial model discussed earlier with 
overshoot stresses felt for the duration of an avalanche, 
once all segments have moved at least once, an increase 
of the force by M/Z is needed to cause another avalanche 
and when one does occur it is likely to run away after its 
diameter becomes of order the velocity correlation length 
£ v of the moving phase at that force. Why this is the con- 
trolling length scale is explained in the next section. 

As discussed in the previous subsection, the pecu- 
liar data for avalanche distributions etc, for f™ ax = 
0.5, M = 0.6 appear to indicate that in the regime 
F > F* the cutoff for the distribution of avalanche areas 
grows as a power of the area. This would suggest the sec- 
ond scenario as such a crossover in the distribution is pre- 
sumably associated with a tendency of larger avalanches 
to runaway. In contrast, the avalanche distributions ob- 
served in this regime with strong pinning suggest that 
there may be a crossover from dissipative- critical-like be- 
havior for smaller sizes to dissipative- cumulative-like be- 
havior for larger sizes, with the crossover occurring for 
y/a of order 0.1 — 0.3L; if this fixed fraction of L persists 
for larger sizes, it would support the third scenario. 

In the next section we develop these ideas further and 
test them in the context of the restarting of the manifold 
after it has come to a stop from the moving phase. 

VII. DYNAMICS OF NUCLEATION 

In the previous section we studied some aspects of how 
the system becomes depinned when it has not previously 
undergone macroscopic motion. We now return to the 
related question of how the manifold restarts after it has 
been stopped by a slow decrease of the force from the 
moving phase. 

As discussed earlier, the stopping process leaves be- 
hind an at-least-somewhat depleted region of small local 
forces. This means that as the force is increased back 
up, the rate of avalanche production, r(F), (which in an 
infinite system is just the probability density of the lo- 
cal forces at zero) will be small. In Fig. |25|(d), r(F) is 
plotted as a function of F for manifolds that have previ- 
ously been moving. Even with relatively strong pinning, 
f™°* = 1.0, it is seen that for M = 0.8 and M = 1.0 
the avalanche production rate is almost three orders of 



magnitude lower than its typical magnitude in a never- 
moved dissipative system. Yet r(F) is relatively flat over 
the regime in which almost all systems will restart: from 
F} to F} + M/Z. Note, however, that in the upper end 
of this range, at least for the data we have taken for 
L = 128, there are no samples still pinned. For smaller 
samples, the range over which the distribution of FJ h st 
extends is broader, i.e. from Fg to Fj + M/Z. Note 
that for all these strong pinning samples, the range over 
which the restarting occurs is two orders of magnitude 
larger than the width of the distribution of Fj . 

A. Distribution of restarting avalanches 

Although some fraction of the samples of size L = 128 
will restart the first time a site is triggered, by no means 
all will and in general there will be a distribution of 
avalanche sizes before restarting. In Fig. |l^(b) the dis- 
tribution of the areas of these avalanches are shown for 
the same parameters as above. The data are somewhat 
sparse due to the small number of avalanches that typi- 
cally occur in a given sample, yet some are observed out 
to about a tenth of the system area. Over the observed 
range, these avalanche distributions are, perhaps surpris- 
ingly, quite close to power laws, with exponents close 
to 0.9 This should be compared with the much smaller 
avalanche area distribution exponent k/2 w 0.39 of the 
critical region avalanches in the dissipative limit. But the 
measured exponent it is close to that of the cumulative 
distribution of dissipative avalanches. 

B. Bubble nucleation 

We conjecture that in the stopped states the length 
scale above which an avalanche is likely to runaway and 
restart the overall motion is of order the velocity corre- 
lation length, £ V (F), of the moving phase at the same 
force. This can be rationalized from the behavior of the 
moving phase. 

At a force that is only a small amount / above F*, 

f = F-F$, (70) 

the moving phase is strongly fluctuating on scales smaller 
than its correlation length, ~ -p. On scales smaller 
than this, there is substantial starting and stopping and 
the motion is fractal. Given the M independence of many 
of the qualitative properties of the moving phase, the mo- 
tion on scales smaller than is probably qualitatively in- 
distinguishable from a finite dissipative avalanche. Thus 
a reasonable guess is that an avalanche from the stopped 
state at / will be fractal on small scales and is not likely 
to runaway unless it becomes of order £„(/). But if it 
does become bigger than this, it will "think" that it is 
part of the moving phase and its local velocity will be un- 
likely to fluctuate to zero. As its local motion continues, 
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it will impart higher stresses to neighboring regions and 
its size will grow. With its interior moving at approxi- 
mately the steady state velocity v(F), we would expect 
the diameter of the resulting "bubble" to grow linearly 
in time. 

In order to test this hypothesis, we investigate in our 
simulations how the manifold begins to move when the 
applied force is increased after it has been stopped. For 
j-max _ 5 and M = o 2j i t can b e seen f rom F ig. || 

that restarting will typically occur at F" c T hyst = F} + M/Z 
by which point the steady state velocity of the moving 
phase is already substantial and its correlation length 
short. In Figs. ^ (a-1) the temporal evolution of the 
restarting process is shown; we see that it is not fractal, 
but looks more like an expanding nucleation bubble as 
occurs after supercooling through an equilibrium first- 
order phase transition. After a short initial transient, 
there is a front that propagates outward from an initial 
seed in a roughly deterministic manner. We observe that 
the radius of the nucleation bubble R(t) grows approxi- 
mately linearly with time so that R(t) — ct, c thus being 
an expansion rate. The interior of the bubble is, as ex- 
pected, moving at roughly the steady state velocity, v(F) 
of the moving phase. [Note that c is a speed in the spatial 
rather than that in the displacement (h) direction which 
is v.] 




FIG. 26. 



Distribution of Fj h t 's with M 



1 0.515 

c hyst 

: 0.2 and L = 64 

comparing f™ ax = 1.5 and f™ ax = 0.5. The rms of F C L for 
the these parameters are 0.4537 ± 0.0009 and -0.139 ± 0.001 
respectively. The up arrows on the horizontal axis indicate 
< Fj > +M/Z, the maximum f™ ax . 



The expansion rate can be estimated from the time 
dependence of the spatially averaged velocity, v{t) by as- 
suming a circular bubble of radius r(t) — ct so that 
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as a function of the time t after the restarting was trig- 
gered, for several different M's and /™ a:E 's. In general, 
at short times after nucleation, fractal growth that is like 
that in a dissipative avalanche occurs within an errati- 
cally growing region of typical diameter £(t) ~ t~ with 
resulting cumulative moment of the avalanche: m(t) ~ 

l d+ t. This results in an apparent c(t) ~ (tkiw)^ ~ 
~ 7 ^ r . See Fig. |^ (b). When the size of the 
avalanche has reached £ V (F), c(t) stops decreasing and 
the flatfish parts of the curves at intermediate times are 
indicative of approximately uniform expansion; the in- 
ferred c(t) in this regime is roughly the expansion rate, 
c, of the bubble. If the expansion rate is small, as will 
occur if the steady state velocity is low, then we expect 
c ~ ^r- ~ £j — * w itb t v the characteristic relaxation time 
in the moving phase; this is also of order the time at 
which the crossover from fractal to bubble growth oc- 
curs. 

At long times, the whole system is moving and c(t) ~ 
as shown. This saturation occurs at t ~ L/c. Be- 
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cause the region near the origin of the bubble has already 
been moving for some time when the regions far away 
start moving, at the time when the furthest regions start 
moving, the roughness of the manifold will be greater 
than it is in steady state The elasticity will diffusively 
smooth out this roughness at longer times leaving be- 
hind only the logarithmic roughness characteristic of the 
moving phase. This excess roughness near the crossover 
time from bubble nucleation to steady state motion is 
the cause of the peak in the width, w 2 (t), of the mani- 
fold shown in Fig. || Naively, we expect the peak w 2 (t) 
to scale as system diameter, L, but subtleties associated 
with elastic slowing down of the velocity inside the bub- 
ble might need to be taken into account to understand 
the process more fully. 

For fp nax — 0.5, the forces at which the restarting oc- 
curs are substantially above Fj and bubble growth ap- 
pears to occur with an expansion rate for M = 0.6 of 
c ~ 0.7, close to the maximum possible rate of unity, 
while for M — 0.2, it is somewhat slower: c w 0.5. In 
the former case, the motion inside of the bubble is essen- 
tially alternating sublattices, behavior characteristic of 
the plateau in the steady state velocity at v « \ that oc- 
curs with substantial overshoots; see Fig. |l|(a). For the 
stronger randomness, f™ ax — 1.0, the hysteresis loops 
are narrower and the nucleation occurs to a lower veloc- 
ity state. For M = 0.2, the expansion rate is quite small, 
c ps 0.13 but the flat constant c regime in the figure is 
observable over a factor of five in time. For these pa- 
rameters, at early stages the evolution appears fractal, 
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consistent with expectations. For M = 0.8, fp iax = 1.0, 
there is a much narrower flat region and the evolution is 
not obviously bubble-like. This can be seen in the snap- 
shots of the local motion as the bubble expands shown in 
Figs. [2^. At long times, all of the samples show, as they 
must, approach of the spatially averaged velocity to its 
steady state value: this asymptotic behavior is indicated 
in Fig. |27|(a) by the solid line. We see that the time at 
which v(t) becomes close to v is roughly that that would 
be expected for a bubble to encompass the whole system 
and run into its periodic "images" . 

More investigation of the dynamics of nuclcation, in 
particular of the transition from fractal on small scales to 
bubble-like on larger, and the dependence of the crossover 
scale on system size, on F — F~, and on the other param- 
eters, is clearly needed. 




1 1 1 1 1 1 1 1 1 1 1 — 

100 200 300 400 500 

t 

(b) 




1 



FIG. 27. (a) Effective bubble growth 

rate, c(t) = ( L n £$ ) 1 ^ 2 as a function of time averaged over 
8 L = 128 restarting samples. The solid line is (^p-) 1//2 , the 
long-time steady state behavior, (b) Log-log plot at early 
times. The solid line is the theoretical expectation in the 
fractal regime with a slope of 4^2. - | w —0.64. 
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t=20, M=0.2, fpm=O.S 



(d) 



t=50, M=0.8, fpm=1.0 
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C. Nucleation avalanche statistics 
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FIG. 28. Dynamics of nucleation: segments that are mov- 
ing at each of the indicated times t after nucleation are shown. 



Armed with the picture of avalanches in the hys- 
teresis loops as fractal on scales out to of order £ V (F) 
and bubble-like on larger scales, we can understand the 
key aspects of the avalanche statistics in this regime. 
Avalanches that are triggered very close to F* — at forces 
within the distribution of F^ for that system size — will 
be fractal on all scales and very unlikely to run away. We 
expect their statistics to be similar to those that occur in 
the regime near Fj when the force is increased from an 
initially flat configuration: the probability of their area 
being larger than a will decay only as a~a out to the 
largest sizes. But as the force is increased, the small 
rate of avalanche production means that there is still a 
substantial probability that no system sized avalanche 
will have occurred until a force well outside the distri- 
bution of F^ . In this regime, avalanches that do occur 
will still have a probability of reaching area a that de- 
cays as out to scales of order £ V (F ), but this is now 
considerably smaller than the system size. Those that 
do reach areas of order £J will tend to runaway and so 
the distribution of areas, p(a, F)da, will be cutoff at this 
independent crossover scale. The probability of runaway 
of an avalanche triggered at F is roughly the probability 
that its area become larger than a v (F) ~ £ v (F) d : 

P Iua (F)^l- \(a,F)da~^-^ (73) 

where p(a, F)da is the fraction of avalanches that occur 
near F which have areas within da of a. In order to 
obtain the cumulative distribution of avalanches in the 
hysteresis loop we must multiply the distribution p(a, F) 
by the probability that runaway has not yet occurred, 

,F 

P P inncd(F, L) = exp ( — / r{F')L d P run (F r )dF') (74) 

and then integrate this over F. If the rate of avalanche 
production per unit area were of order unity, then P p i nn ed 
would be small for F — F} > L~i and the cumulative 
avalanche distribution would be dominated by avalanches 
in the regime very near Fjr for which the distribution is 
characteristic of the critical point; the cumulative distri- 
bution of avalanche areas in the hysteresis loop would 
then decay with the exponent n/d consistent with what 
we observe for the dissipative case. But with r(F) very 
small for strong pinning and M = 0.8 or 7\f = 1.0 — al- 
most three orders of magnitude smaller than for the dis- 
sipative case, (see Fig. ||](d)) P P mned will only become 
small when < rL d corresponding to F—F^ > r~ 
a factor of fifty or so further from F} than in the dissi- 
pative case, roughly consistent with what is observed. 

VIII. CROSSOVERS AND ASYMPTOPIA 
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A. Hysteresis loops 

The form of the avalanche distribution in hysteresis 
loops with a very low density of nucleation segments 
should be apparent from the discussion at the end of the 
previous section: Most avalanches will take place when F 
is close to the -Fj at which the system restarts, thus it 
is the distribution p(a, (Fj h ))da in this regime that will 
dominate the cumulative statistics. This results in a dis- 
tribution of avalanche areas that for small avalanches is 
critical with exponent § , but for large avalanches decays 
with the larger cumulative exponent K cum = 1. With 
small r, as occurs in the strong pinning systems with 
M = 0.8 and M — 1.0, the crossover in a system of 
linear size L occurs at a length scale of order Lr~3 or 
corresponding crossover area of order rL 2 . For the size 
ranges investigated here — up to L = 128 — this means 
that in practice almost the full range of avalanche areas 
will be in the large avalanche tail region with exponent 
J( cum = 1 . But extrapolating to very large systems from 
these data would be highly misleading: in the limit of 
large systems, the distribution should scale like the crit- 
ical dissipative limit over most of the range of avalanche 
sizes with the rapid decay observed in our data obtaining 
only for the largest few orders of magnitude of avalanche 
areas. 

The dynamics of a typical runaway event will reflect 
this crossover. Initially, it will start as a fractal avalanche 
similar to those in the dissipative limit. But when — and 
if — an avalanche reaches the crossover length of order 
Lr d , it will crossover to bubble like expansion with its 
interior behaving like the moving phase and its radius 
expanding linearly in time at a rate, c, that scales as 
(Lrd) 1 ^ 2 . As the distribution of fj hyst — F} is broad, 
however, there will be substantial run-to-run variations 
in both the crossover size and the rate of growth of the 
bubble. 

The behavior of the size dependence of the distribu- 
tion of widths of the hysteresis loops, as shown in Figs. 
|TT|(a) (b), can also be better understood in this frame- 
work. In small systems, the force will typically have to be 
increased by M/Z from F* in order to trigger restarting. 
In contrast, in very large samples, we expect that the 
restarting will typically occur at forces F = (Fj) + gL~~ 
which are above (F}) by of order the width of the dis- 
tribution of F*. Yet because of the strong suppression 
of the density, r(F), of seeds for nucleation of avalanches 
in the stopped system, in practice the numerical factors, 
g can be very large. The distribution of g in the small 
r limit is broad, and has a shape close to a Weibull dis- 
tribution with characteristic scale l/r~ and shape pa- 
rameter 4j but we do not expect it to be exactly of this 
form. The reason is due to the breakdown of the argu- 
ment given in Section V: In sufficiently small subregions 
of size (-jj 1 ) of a system of size L d , the possible nucleation 
processes are not independent from one subregion to the 
next. This is because the diameter of a critical nucleation 



avalanche when it crosses over to approximately deter- 
ministic bubble- like growth is of order £ v (F) which is of 
order a fixed — albeit r dependent — fraction of L at the 
typical forces at which the first runaway avalanche occurs 
in a large system. Thus the division into approximately 
independent subregions will not work for b > 1/r^. 

Because the hysteresis loops vanish in the limit of large 
systems, the total number of avalanches per unit area 
that will occur in a hysteresis loop decreases as a power of 
system size. This will also occur in the dissipative limit, 
however in that case the "runaway" avalanches that occur 
for F rs F^ will not have all that much larger moment 
than ones that only involve, say, half the system area. 
But due to the strong suppression of avalanche produc- 
tion after a large event, in the presence of substantial 
stress overshoots the total number of avalanches will be 
further reduced in spite of the hysteresis loop becoming 
wider. The above analysis implies that the total number 
of avalanches per unit area that occur in the hysteresis 
loop will be of order: 

n a val ~ r^^-^j- . (75) 

L u 

Furthermore, the runaway avalanche will typically run 
for a very long way — probably exponentially far in 1/r 
- before stopping at a displacement, h, for which the 
total pinning is anomalously strong. 

At this point, it is not clear that the numerical data are 
consistent with this picture. In particular, it appears that 
the size dependence of the widths of the hysteresis loops 
is somewhat weaker than would be expected asymptoti- 
cally. Nevertheless, we believe that the deviation of the 
scaling of the width of the hysteresis loops from the L~~ 
scaling expected asymptotically is most likely associated 
with the crossover from the small system to large system 
behavior. In sufficiently large systems, we conjecture that 

H = v . (76) 



B. Initial depinning 

The complicated avalanche behavior observed when 
depinning is approached for the first time starting from 
almost flat initial conditions is probably due to a hy- 
brid of effects. For the strong pinning samples with 
substantial-sized overshoots, in systems of size 128 2 
enough avalanches typically occur before runaway that a 
strong suppression of the avalanche production rate oc- 
curs over the whole system. This suppression enables the 
force to be increased well past the range in which run- 
away would occur in the dissipative limit, although not 
as far as it can for samples that have been stopped from 
the moving phase. The statistics of the avalanches in 
this regime will have a similar form to that of avalanches 
in the hysteresis loops discussed above, with a crossover 
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from critical-like to cumulative-like at a length scale that 
depends on how strongly r(F) has been suppressed. 

For the weaker pinning samples, for example with 
frnax _ q 5 anc j m _ q.6, the rate of avalanche produc- 
tion is not suppressed nearly as strongly and it is most 
likely that the crossover regime is right in the middle of 
the range of sample sizes studied. 

Further systematic tests to investigate this overall sce- 
nario would be very useful. In particular, one would like 
to understand on what aspects of the model, parameters, 
and history, crucial quantities such as the suppression of 
avalanche production rate after a large avalanche, de- 
pend. 

IX. DISCUSSION AND CONCLUSIONS 

In this paper we have introduced and studied numeri- 
cally and by scaling analyses a simple model of dcpinning 
transitions in which there are dynamic stress overshoots 
caused by rapid motion of segments of the manifold. 

A. Macroscopic behavior 

We have produced a substantial amount of evidence 
that the macroscopic behavior of this model with var- 
ious kinds of generic initial conditions is not hysteretic 
for a wide range of stress overshoots in spite of the 
coexistence over a range of the driving force of static 
configurations and steady state moving configurations. 
In particular, the critical force, F C (M) is uniquely de- 
fined and history independent in the limit of large sys- 
tems; the distributions of the critical force, F}, at which 
steady state motion ceases, that at which it restarts when 
the force is increased from such a stopped state, -Fj hyst , 
and that at which macroscopic motion commences when 
the force is increased from a flat initial condition for 
the first time, Fj , are all narrow and converge to 
F C (M) in the limit of large systems. Furthermore, the 
steady state velocity, v(F), is a unique function of the 
driving force, F that vanishes as a universal power of 
F — F c . Nevertheless, in principle there arc static config- 
urations that could be reached by careful enough control 
of the dynamics of all the system for all forces up to 
F cQ = F C (M = 0) > F C (M). But if the system were to 
be started in such a configuration and the force increased 
by an amount that is arbitrarily small in the limit of a 
large system, an avalanche would be triggered that would 
runaway and cause the whole system to start moving. 

The universality class of the dcpinning transition for 
all M less than a multicritical value M c appears to be 
that of the purely dissipative limit without stress over- 
shoots that has been analyzed by renormalization group 
methods previously. This universality class is thus much 
broader than had been conjectured. 



The crucial feature of the dynamics on which this 
macroscopic uniqueness and universality relies is the ex- 
istence of a non-zero density of nucleation sites for new 
avalanches after an avalanche has run through a region. 
These must exist for arbitrarily small increases of the ap- 
plied force even if their density, r, is much lower than it 
would have been in the absence of stress overshoots. If 
these did not exist — more precisely if there were a deple- 
tion layer of local forces after an avalanche — then there 
would be macroscopic hysteresis with the force needing to 
be increased by the width of the depletion layer to restart 
the system once it has stopped after having moved pre- 
viously. 

B. Finite-size effects 

In the remainder of this last section we summarize the 
behavior that occurs in our model for a range of M in 
which there is a strongly suppressed density of nucleation 
seeds. We also consider briefly under what circumstances 
similar behavior will obtain; and discuss some of the con- 
sequences for a system in which understanding the finite 
size effects is crucial: earthquakes on a geological fault. 

In the dissipative limit of our model, all finite size 
effects in the vicinity of the depinning transition occur 
within a range of forces around F c of width of order 1/L~ . 
In contrast, in the presence of stress overshoots finite- 
size effects can be important over a much wider range 
of forces. These effects will be particularly pronounced if 
the suppression of nucleation seeds for further avalanches 
is very strong, i.e. if after a system spanning avalanche 
the density of nucleation segments, r, is very small. 

With r small, unless samples are sufficiently large, 
L d > Z/(Mr) there will usually be no nucleation un- 
til the force has been increased by the magnitude, M/Z, 
of the overshoots. Unless M is very small, the nucle- 
ation will then be bubble- like with rapid motion in a lin- 
early expanding bubble occurring as soon as one segment 
is triggered. For larger systems, the width of hysteresis 
loops will start to decrease and eventually, if the scenario 
presented above is correct, decrease as 1/L~ but with a 
large coefficient that is a randomly varying multiple of 
l/r~ . Sites that are triggered as F is increased after 
the whole system has been stopped may result in finite 
avalanches with a distribution that crosses over from one 
power law to another as a function of their size. Eventu- 
ally one of these avalanches will runaway: initially it will 
grow in a fractal manner, but when its diameter reaches 
of order L/r~d } it will began to expand in a more deter- 
ministic manner at an approximately constant rate until 
it covers the whole system and the macroscopic motion 
restarts. In the limit of small r this regime of bubble 
growth will obtain over a wide range of length scales. 
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C. Model earthquake dynamics and statistics 

The picture presented here for effects of stress over- 
shoots on the dynamics of finite size systems that occur 
as the driving force is varied have particularly interest- 
ing consequences for models of earthquakes on disordered 
faults. The appropriate driving to model a geological 
fault is not a constant force, by rather driving by a weak 
spring. This is roughly equivalent for our case with short 
range elastic interactions to replacing the applied force 
by a "pulling spring" with 
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F(t) = j^{v s t-h(t)) 



(77) 



where G is an effective clastic constant of order unity, the 
1/L 2 factor arises from the V 2 elasticity, and 



h(t) 



(78) 



is the spatially averaged displacement. 

As has been discussed previously, in the absence 
of stress overshoots a system driven with an infinitesi- 
mal "shearing" velocity, v s , "self-organizes" into a sta- 
tistically steady state with a power law distribution of 
avalanche diameters that falls off as -p^j- out to length 
scales of order L, (the cutoff being affected by the magni- 
tude of G) . There is no particular qualitative distinction 
between the character of "earthquakes" that have sub- 
stantial I -C L and those with I ~ L. 

With stress overshoots that give rise to a substan- 
tial reduction of the density of seeds for nucleating new 
events, the behavior is quite different. As the pulling 
spring advances there will be a string of avalanches until 
one runs away and the whole system moves. This mo- 
tion will only stop when the force F(t) has decreased as 
a result of the motion to , a force that will have varia- 
tions of order 1/L~ . As the spring slowly restretches, the 
pulling force will increase again and a series of quakes of 
various sizes will occur at a rate given by rL d /v s . These 
events will not relieve much of the accumulating stress 
and so the force will continue to increase until a force 
f J h st at which a runaway event occurs. The typical mag- 
nitude of FX — F} and the width of its distribution are 

<--hyst *- 

both of order -V^j, The evolution of the runaway event 

will at first be fractal, but then a crossover to bubble-like 
growth will occur and it will grow with constant expan- 
sion rate until the whole system is moving. 

The initial slip velocity of the large earthquake will be 
of order v(Fj ) which will vary substantially, but have 
typical magnitude 
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As it runs, the driving force will gradually decrease until 
it gets down to another F^ at which point it will stop. 



The total displacement in such an event will be of order 
Ah ~ L 2 (Fj h t — F}) and its moment hence of order 



m R ~ L d Ah 
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bigger by the r dependent factor than the typical largest 
events in the dissipative limit. [In obtaining this result, 
the scaling law £ + i = 2 was used.] 

The distribution of earthquake moments in this model 
will be a composite of characteristic earthquake and 
Gutenberg-Richter-like statistics p2[ , |23|. Small events 
will have a power law distribution with exponent n/(d+Q 
that will obtain out to moments of order 
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much smaller than m,R for small r, The distribution of 
larger events will fall off more rapidly with exponent 
d/ {d+C) out to events with moments of order ra ~ L d+l > . 
Events larger than this will not occur unless they are run- 
away events. As these run for a considerably larger Ah 
before stopping, there will be a strong suppression of the 
distribution of moments from m Q out to of order uir, but 
the runaway events will have a non-trivial distribution. 

This picture implies interesting anticorrelations be- 
tween the occurrences of the largest of the intermediate 
size and the runaway events: the former are more likely 
to occur relatively soon after a runaway event rather than 
just before one because at later stages in the "earthquake 
cycle" , once an event gets to a small fraction of the sys- 
tem size, it will almost certainly runaway. Some infor- 
mation on the statistics of the largest events can also 
be inferred from the approximate Weibull form of the 
distribution of FJ h st which will primarily determine the 
moments of the largest events. 

In order to translate these results to models more ap- 
propriate for geological faults, the elasticity needs to be 
long-ranged, with static stress transfer falling off as 



■/(x,y) 



l x -yr 



(82) 



in the two dimensional case of interest for all but the 
biggest earthquakes. This modifies scaling laws to £ + 
i = 1 and k = d — 1 + (. In the absence of long range 
correlations in the random properties of the fault, in the 
dissipative limit the predicted value of £ is zero, up to 
logarithms. But with long-range correlated randomness 
it could be substantially larger. 

The consequences of the scenario suggested by our re- 
sults for geological fault dynamics certainly merits fur- 
ther investigation along these and other lines. But to 
conclude this paper, we turn to another important issue: 
how much can be carried over from the stress overshoot 
model to other types of local dynamic effects? We should 
remark that this scenario is valid within the context of 
fi = v for very large systems. 
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D. Other dynamic overshoot effects 

At the beginning of this paper we argued that some 
of the effects of local inertia of a driven, pinned elastic 
manifold would be qualitatively similar to the effects of 
stress overshoots. A key question to raise at this point 
is whether the finite density of seeds for future avalanche 
events will be left behind in inertial systems. We con- 
jecture that this will be the case if the inertia is not too 
large, as we found here for stress overshoots, but this 
bears more careful thought. 

In the context of earthquakes, a more important is- 
sue is dynamic frictional weakening. An extreme form 
of frictional weakening was considered in |L6| |: if once a 
segment has moved it is always easier for it to move again 
at any time in the same earthquake, this will certainly 
give rise to macroscopic hysteresis. As long as there is 
healing back up to a higher strength between events, a 
system spanning event will leave behind a configuration 
that cannot start moving again until the driving force 
has increased again by enough to overcome the difference 
between the static and dynamic friction. In reality, the 
history dependence of frictional forces is complicated and 
some healing will occur already during an event on those 
segments that slow down or stop while other parts are 
still moving. But whether this, combined with dynamic 
stress overshoots, will leave behind a finite density of nu- 
cleation sites for easily triggerable events is a difficult but 
important question. 

Although this paper has perhaps raised as many ques- 
tions as it has answered, the progress it represents in un- 
derstanding of the dynamics of driven elastic manifolds 
in the presence of both randomness and non-dissipative 
dynamics should help frame and address some of the key 
questions that remain. 
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